Простейшие модели экономической динамики

604 days ago by kpp

def myint(f, start, h, N): yield start for _ in range(1, N): start = start + h*f(start) yield start 
       
@interact def lfi(s=2+2/10, a=-9/10, T=4, N=50, interval = input_box([0, 2]), ): y = var('y') S = s*(y+a)**3 - a**3 + 2 I = 3*y + 5/2 - S f = I - S dfdy = f.differentiate() html('<h3>Логистические функции потребления и инвестиций</h3>') html('<center><font color="blue">$S(y) = %s$</font></center>'%latex(S)) html('<center><font color="green">$I(y) = %s$</font></center>'%latex(I)) p1 = plot(S, interval, color='blue', axes_labels=['$y$', '']) p2 = plot(I, interval, color='green') show(p1 + p2) html('<center><font color="red">$f(y) = %s$</font></center>'%latex(f)) html(r'''<center><font color="black">$\frac{df}{dy} = %s$</font></center>'''%latex(dfdy)) p3 = plot(f, interval, color='red', detect_poles='show', axes_labels=['$y$', '']) p4 = plot(dfdy, interval, color='black') show(p3+p4, ymin=-3, ymax=3) assume(y, 'rational') yroots = f.solve(y, to_poly_solve='force') yroots = [y.rhs().n(digits=16) for y in yroots] yroots.sort() unstable = None; s = r'''<table width="100%"><th> <tr><td>$y_{root}$</td> <td>$S(y_{root})$</td> <td>$\frac{df}{dy}$</td> <td>устойчивость</td> </tr></th>''' for yroot in yroots: s = s + "<tr>" point = yroot#.rhs().n(digits=5) val = S(y=point) dval = dfdy(y=point) u = "не устойчивая"; if dval < 0: u = "устойчивая" else: unstable = yroot tab = ["<td>$"+str(i)+"$</td>" for i in [point, val, dval, u] ] for t in tab: s = s + t s = s + "</tr>" s = s + "<table>" html(s) h = T/N; fnum = lambda t: n(f(y=t)) points = [k*h for k in range(0,N)] starts = [unstable + 0.5, unstable - 0.5, unstable + 0.001, unstable - 0.002] plots = [] for i in range(0, len(starts)): clr = hue(float(i)/len(starts)) pt = starts[i] plots.append(list_plot(zip(points, myint(fnum, pt, h, N)), color=clr)) p8 = plot(yroots, 0, T, linestyle='--', color='black') #show(p5+p6+p7+p8) show(sum(plots) + p8) 
       

Click to the left again to hide and once more to show the dynamic interactive window

html("<h2>Динамика несвязанных секторов экономики</h2>") @interact def _(a0 = input_box( [1+2/1000, 98/100, 1], label="$a_0 = $"), a1 = input_box( [5382/10000, 6/10, 4/10], label="$a_1 = $"), a2 = input_box( [4518/10000, 4, 5], label="$a_2 = $"), alpha = input_box( [2/10, 3/10, 4/10], label=r'$\alpha = $'), beta = input_box( [1/10, 2/10, 3/10], label=r'$\beta = $'), h = 1, N = 100, K0 = input_box([1,1,1], label="Начальная капитализация: "), L0 = input_box([1,1,1], label="Трудовые затраты: "), ): html("<h3>Динамика изменения капитала по секторам</h3>") K0 = vector(K0) L0 = vector(L0) K = var('K') YY = lambda K,L,i: a0[i]*K[i]**a1[i]*L[i]**a2[i] Y = lambda K,L: vector([YY(K,L,i) for i in range(0, len(K))]) Y = fast_float(Y) A = diagonal_matrix(alpha) B = diagonal_matrix(beta) F = lambda K: (A*Y(K,L0)-B*K).n() points = [k*h for k in range(0,N)] values = [value for value in myint(F, K0, h, N)] plots = [] for i in range(0, len(K0)): c = hue(float(i)/len(K0)) plots.append(list_plot(zip(points, [v[i] for v in values]), color=c ) ) show(sum(plots)) vechtml = lambda vec, name: html( "$$%(var)s = %(val)s$$"%{'var':name, 'val':latex(vec.transpose())} ) html("Величина капитала в конце периода по секторам:") KN = values[-1] vechtml(KN, 'K_{fin}') html("Величина дохода в конце периода по секторам:") YN = Y(KN, L0) vechtml(YN, 'Y_{fin}') html("Суммарный доход: $%f$."%sum(YN)) al = diagonal_matrix(alpha) I = al*YN html("<br />Расходы на капитализацию по секторам:") vechtml(I, 'I') html("Расходы на капитализацию в целом: $%f$."%sum(I)) S = YN - I html("<br />Потребление по секторам:") vechtml(S, 'S') html("Потребление в целом: $%f$."%sum(S)) 
       

Click to the left again to hide and once more to show the dynamic interactive window

html("<h2>Динамика связанных секторов экономики</h2>") @interact def _(a0 = input_box( [1+2/1000, 98/100, 1], label="$a_0 = $"), a1 = input_box( [5382/10000, 6/10, 4/10], label="$a_1 = $"), a2 = input_box( [4518/10000, 4, 5], label="$a_2 = $"), A = input_box( [[1/10,0,0], [1/10,2/10,0], [0,1/10,0]], label="$A = $"), alpha = input_box( [2/10, 3/10, 4/10], label=r'$\alpha = $'), beta = input_box( [1/10, 2/10, 3/10], label=r'$\beta = $'), h = 1, N = 100, K0 = input_box([1,1,1], label="Начальная капитализация: "), L0 = input_box([1,1,1], label="Трудовые затраты: "), ): K0 = vector(K0) L0 = vector(L0) A = matrix(A) html("$$A = %s$$"%latex(A)) html("<h3>Динамика изменения капитала по секторам</h3>") K = var('K') YY = lambda K,L,i: a0[i]*K[i]**a1[i]*L[i]**a2[i] Y = lambda K,L: vector([YY(K,L,i) for i in range(0, len(K))]) Y = fast_float(Y) B = diagonal_matrix(beta) F = lambda K: (A*Y(K,L0)-B*K).n() points = [k*h for k in range(0,N)] values = [value for value in myint(F, K0, h, N)] plots = [] for i in range(0, len(K0)): c = hue(float(i)/len(K0)) plots.append(list_plot(zip(points, [v[i] for v in values]), color=c ) ) show(sum(plots)) vechtml = lambda vec, name: html( "$$%(var)s = %(val)s$$"%{'var':name, 'val':latex(vec.transpose())} ) html("Величина капитала в конце периода по секторам:") KN = values[-1] vechtml(KN, 'K_{fin}') html("Величина дохода в конце периода по секторам:") YN = Y(KN, L0) vechtml(YN, 'Y_{fin}') html("Суммарный доход: $%f$."%sum(YN)) al = diagonal_matrix(alpha) I = al*YN html("<br />Расходы на капитализацию по секторам:") vechtml(I, 'I') html("Расходы на капитализацию в целом: $%f$."%sum(I)) S = YN - I html("<br />Потребление по секторам:") vechtml(S, 'S') html("Потребление в целом: $%f$."%sum(S)) 
       

Click to the left again to hide and once more to show the dynamic interactive window

@interact def _(a0 = input_box(1+2/1000, label=r'$a_0 = $'), a1 = input_box(5382/10000, label=r'$a_1 = $'), a2 = input_box(4618/1000, label=r'$a_2 = $'), alpha = input_box(2/10, label=r'$\alpha = $'), beta = input_box(1/10, label=r'$\beta = $'), K0 = input_box(10, label=r'$K_0 = $'), eps = input_box(1/100000, label=r'$\varepsilon = $'), ): K,L = var('K,L') Y = a0 * K**a1 * L**a2; html(r'$$Y = %s$$'%latex(Y)) f = alpha * Y - beta * K html(r'$$f = %s$$'%latex(f)) DfDK = f.differentiate(K) html(r'$$\frac{\partial f}{\partial K} = %s$$'%latex(DfDK)) eq = (f(L=1) == 0) html('$$%s$$'%latex(eq)) assume(K,'rational') assume(K>0) c1 = f(L=1)._maxima_() c1.eval("load(newton1);") Kroot = float( c1.newton(K, K0, eps) ) html('$$K_{root} = %s$$'%latex(Kroot)) html(r'$$\left.\frac{\partial f}{\partial K} \right|_{K_{root}} = %s$$'%\ latex(DfDK(K=Kroot, L=1))) 
       

Click to the left again to hide and once more to show the dynamic interactive window