CompEc 2022/2023

Computational Economics
Università degli Studi di Firenze

SYLLABUS

Official Syllabus from University of Florence website

Vai al link

A TENTATIVE PROGRAM

Vai al link

1. Income Expenditure Model

Notes
Python Scripts:

Income-Expenditure Model (Graph)

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Sun Mar 5 2023 @author: Giorgio Ricchiuti @website: www.grarchive.net Description Income-Expenditure Model """ import numpy as np import matplotlib.pyplot as plt """ Parameters """ C1=10000 #autonomous consumption c=0.8 #marginal propensity to consume I=2500 # investments G=3000 #public expenditure X=500 #exports Q1=200 #autonomous imports q=0.15 #marginal propensity to import t=0.15 #tax rate TA1=100 #fixed taxes TR=50 #transfers """ Auxiliary variables """ A=C1+c*(TR-TA1)+I+G+X-Q1 # sum of the autonous components a=c*(1-t)-q # slope of the expenditure line """ Graph """ plt.title("Income-Expenditure Model") Y=np.linspace(0,45000,101) Z=A+a*Y plt.grid(True) plt.plot(Y,Z) plt.plot(Y,Y) # an increase in autonomous consumption C1=5000 # Definition/Function to produce a plot def plot_IE(A,a): Y=np.linspace(0,45000,101) Z=A+a*Y plt.grid(True) plt.plot(Y,Z) plt.plot(Y,Y) plot_IE(A,a)

Static Analysis

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Sun Mar 5 2023 @author: Giorgio Ricchiuti @website: www.grarchive.net Description Income-Expenditure Model: static analysis """ import sympy as sym """ Parameters """ C1=10000 #autonomous consumption c=0.8 #marginal propensity to consume I=2500 # investments G=3000 #public expenditure X=500 #exports Q1=200 #autonomous imports q=0.15 #marginal propensity to import t=0.15 #tax rate TA1=100 #fixed taxes TR=50 #transfers """ Auxiliary variables """ A=C1+c*(TR-TA1)+I+G+X-Q1 # sum of the autonous components a=c*(1-t)-q # slope of the expenditure line def static_IE(A,a): y=sym.Symbol('y') eqn=sym.solveset(A+a*y-y, y) print(eqn) static_IE(A,a) Check=A/(1-a) print('Check: {}'.format(Check))

Dynamic Analysis

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Sun Mar 5 2023 @author: Giorgio Ricchiuti @website: www.grarchive.net Description Income-Expenditure Model: dynamic analysis """ import numpy as np import matplotlib.pyplot as plt def IE(N, x0, A, a, phi): #define a vector of zeros x=np.zeros(N); #initial condition x[0]=x0; #draw the expenditure and 45 degree line y=np.linspace(0,50000,11) z=A+a*y plt.plot(y,z,color="black", label="z=A+a*y") plt.plot(y,y) plt.grid(True) plt.axhline(color='black', lw=0.5) plt.axvline(color='black', lw=0.5) # for t in range(N-1): x[t+1]=x[t]+phi*(A+a*x[t]-x[t]) plt.plot([x[t],x[t]],[x[t],x[t+1]], color='blue') plt.plot([x[t],x[t+1]],[x[t+1],x[t+1]], color='red') plt.plot([0,max(x)], [0, max(x)], color='black', linestyle='dashed', linewidth=0.5) plt.xlabel('Y') plt.ylabel('Z') plt.title('Dynamic Income Expenditure Model ') plt.show() return x par = {'a' : 0.53, 'A' : 15000, 'x0': 10000, 'N':20, 'phi':3.5} x = IE(**par) plt.plot(x) plt.xlabel('Time') plt.ylabel('Income') plt.title('Time series for Income: convergence to the equilibrium') plt.show()

2. Logistic model

Simple mathematical models with very complicated dynamics (Robert May, Nature 1976)
Python Scripts:

Logistic Function

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Fri Mar 10 11:33:34 2023 @author: giorgio """ import numpy as np import matplotlib.pyplot as plt def logistic_map(N, x0, lamb): x=np.zeros(N); x[0]=x0; #let's make a draw of the logistic function z=np.linspace(0,1,101) y=lamb*z*(1-z) plt.plot(z,y,color="black", label="y=lamb*z*(1-z)" ) #a loop to calculate the evolution for t in range(N-1): x[t+1]=lamb*x[t]*(1-x[t]) plt.plot([x[t],x[t]],[x[t], x[t+1]], color='blue') plt.plot([x[t],x[t+1]],[x[t+1], x[t+1]], color='blue') plt.plot([0,1.2],[0,1.2], color='black', linestyle='dashed', linewidth=0.5) return x ''' Figures with different initial conditions ''' f= logistic_map(35, 0.3, 3.9) plt.title('Logistic Map 1') plt.grid(True) plt.show() g= logistic_map(35, 0.3001, 3.9) plt.title('Logistic Map 2') plt.grid(True) plt.show() plt.plot(f, color='blue') plt.plot(g, color='red') plt.xlabel('time') plt.ylabel('x') plt.title('Time series of x') plt.show() ''' Figures for different values of lambda ''' f= logistic_map(35, 0.3, 1) plt.title('Logistic Map 1') plt.grid(True) plt.show() g= logistic_map(35, 0.3, 3.2) plt.title('Logistic Map 2') plt.grid(True) plt.show() h= logistic_map(35, 0.3, 3.9) plt.title('Logistic Map 3') plt.grid(True) plt.show() plt.plot(f, color='blue', label='1') plt.plot(g, color='red', label='3.2') plt.plot(h, color='green', label='3.9') plt.legend(loc=3) plt.xlabel('time') plt.ylabel('x') plt.title('Time series of x') plt.show() """ All figures with subplots """ fig=plt.figure() ax1=fig.add_subplot(321) x = logistic_map(35, 0.3, 2) ax2=fig.add_subplot(322) plt.plot(x) ax3=fig.add_subplot(323) x1 = logistic_map(35, 0.3, 3.3) ax4=fig.add_subplot(324) plt.plot(x1) ax5=fig.add_subplot(325) x2 = logistic_map(35, 0.3, 3.9) ax6=fig.add_subplot(326) plt.plot(x2) #plt.title('Time series 3') plt.show()

Bifurcation py 1

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio """ import numpy as np import matplotlib.pyplot as plt def logistic_eq(r,x): return r*x*(1-x) # Create the bifurcation diagram def bifurcation_diagram(seed, n_skip, n_iter, step=0.0001, r_min=0): print("Starting with x0 seed {0}, skip plotting first {1} iterations, then plot next {2} iterations.".format(seed, n_skip, n_iter)); # Array of r values, the x axis of the bifurcation plot R = [] # Array of x_t values, the y axis of the bifurcation plot X = [] # Create the r values to loop. For each r value we will plot n_iter points r_range = np.linspace(r_min, 4, int(1/step)) for r in r_range: x = seed; # For each r, iterate the logistic function and collect datapoint if n_skip iterations have occurred for i in range(n_iter+n_skip+1): if i >= n_skip: R.append(r) X.append(x) x = logistic_eq(r,x); # Plot the data plt.plot(R, X, ls='', marker=',') plt.ylim(0, 1) plt.xlim(r_min, 4) plt.xlabel('r') plt.ylabel('X') plt.show() return R, X, x R, X, x=bifurcation_diagram(0.2, 100, 5) R, X, x= bifurcation_diagram(0.2, 100, 10) R, X, x= bifurcation_diagram(0.2, 100, 10, r_min=2.8)

Bifurcation py 2

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio """ import numpy as np import matplotlib.pyplot as plt # Create the bifurcation diagram def bifurcation_diagram(N): step=0.001 x=np.zeros(N) X=np.zeros((int(1/step+1),N)) r_min=0 r_max=3.95 x[0]=0.01 r_range = np.linspace(r_min, r_max, int(1/step+1)) flag=-1 for i in r_range: flag = flag + 1 for t in range(N-1): x[t+1]=i*x[t]*(1-x[t]) #print(x) X[flag:] = x #print(X) T_last = 350 Nr =int(N - T_last) plt.show() for i in range(len(r_range)): plt.plot(r_range[i]*np.ones((T_last)), X[i, Nr:N], linestyle='dashdot', linewidth=0.01, markersize=0.1, marker=',', color='black') return x x=bifurcation_diagram(400)
Bibliography:

Hommes, C. (2013), Behavioural Rationality and Heterogenous Expectations, ed. Cambridge University Press (Chapter 2.)

3. Dynamic ISLM

Notes
Python Scripts:

Dynamic ISLM case a)

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio ISLM Model case a) """ import numpy as np import matplotlib.pyplot as plt def ISLMa(N, Y0, r0, a, Z, b, k, h, M, P, phi): #define a vector of zeros Y=np.zeros(N) r=np.zeros(N) #initial condition r[0]=r0 Y[0]=Y0 z=1/(1-a) #draw IS and LM schedule y=np.linspace(20000,30000,11) r_IS=Z/b-y/(z*b) r_LM=k*y/h-(M/P)/h plt.plot(y,r_IS,color="red", label="IS") plt.plot(y,r_LM,color="red", label="LM") plt.grid(True) # plt.axhline(color='black', lw=0.5) # plt.axvline(color='black', lw=0.5) # for t in range(N-1): Y[t+1]=Y[t]+phi*(Z+a*Y[t]-b*r[t]-Y[t]) r[t+1]=k*Y[t+1]/h-(M/P)/h #plt.plot([0,max(x)], [0, max(x)], color='black', linestyle='dashed', linewidth=0.5) plt.xlabel('Y') plt.ylabel('r') plt.title('Dynamic ISLM Model A') plt.scatter(Y,r, c=np.arange(N), cmap='viridis', zorder=3, s=15) plt.plot(Y,r,color="black", zorder=3, linewidth=0.9) plt.show() return Y, r par = {'a' : 0.53, 'Z' : 15000, 'Y0': 25000, 'r0': 25, 'N': 15, 'phi':0.95, 'k':0.95, 'b':100, 'h':100, 'M':20000, 'P':1 } Y,r = ISLMa(**par) plt.plot(Y) plt.xlabel('time') plt.ylabel('Y') plt.title('Time series for Income: convergence to the equilibrium') plt.show() plt.plot(r) plt.xlabel('time') plt.ylabel('r') plt.title('Time series for Interest rate: convergence to the equilibrium') plt.show()

Dynamic ISLM case b)

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio ISLM Model case b) """ import numpy as np import matplotlib.pyplot as plt def ISLMb(N, Y0, r0, a, Z, b, k, h, M, P, phi, psi): #define a vector of zeros Y=np.zeros(N) r=np.zeros(N) #initial condition r[0]=r0 Y[0]=Y0 z=1/(1-a) #draw IS and LM schedule y=np.linspace(0,100,11) r_IS=Z/b-y/(z*b) r_LM=k*y/h-(M/P)/h plt.plot(y,r_IS,color="red", label="IS") plt.plot(y,r_LM,color="blue", label="LM") plt.legend() plt.grid(True) # plt.axhline(color='black', lw=0.5) # plt.axvline(color='black', lw=0.5) for t in range(N-1): Y[t+1]=Y[t]+phi*(Z+a*Y[t]-b*r[t]-Y[t]) r[t+1]=r[t]+psi*(k*Y[t]-h*r[t]-(M/P)) plt.xlabel('Y') plt.ylabel('r') plt.title('Dynamic ISLM Model B') plt.scatter(Y,r, c=np.arange(N), cmap='viridis', zorder=3, s=15) plt.plot(Y,r,color="black", zorder=3, linewidth=0.9) plt.show() return Y, r par = {'a' : 0.5625, 'Z' : 50, 'Y0': 42, 'psi': 0.4, 'r0': 12, 'N': 150, 'phi':0.2, 'k':0.25, 'b':1.525, 'h':0.5, 'M':8, 'P':1 } Y,r = ISLMb(**par) plt.plot(Y) plt.xlabel('time') plt.ylabel('Y') plt.grid(True) plt.title('Time series for Income: convergence to the equilibrium') plt.show() plt.plot(r) plt.xlabel('time') plt.ylabel('r') plt.grid(True) plt.title('Time series for Income: convergence to the equilibrium') plt.show()

...with Fiscal Policy

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio ISLM Model case b) with fiscal policy """ import numpy as np import matplotlib.pyplot as plt def ISLMb(N, Y0, r0, a, Z, b, k, h, M, P, phi, psi): #define a vector of zeros Y=np.zeros(N) r=np.zeros(N) #initial condition r[0]=r0 Y[0]=Y0 z=1/(1-a) #draw IS and LM schedule y=np.linspace(0,100,11) r_IS=Z/b-y/(z*b) r_LM=k*y/h-(M/P)/h plt.plot(y,r_IS,color="blue", label="IS0") Z=75 r_IS2=Z/b-y/(z*b) plt.plot(y,r_IS2,color="blue", linestyle='dashed', label="IS1") plt.plot(y,r_LM,color="red", label="LM") plt.legend() plt.grid(True) plt.axhline(color='black', lw=0.5) plt.axvline(color='black', lw=0.5) for t in range(N-1): Y[t+1]=Y[t]+phi*(Z+a*Y[t]-b*r[t]-Y[t]) r[t+1]=r[t]+psi*(k*Y[t]-h*r[t]-(M/P)) # plt.plot([0,max(Y)], [0, max(r)], color='black', linestyle='dashed', linewidth=0.5) plt.xlabel('Y') plt.ylabel('r') plt.title('Dynamic ISLM Model B') plt.scatter(Y,r, c=np.arange(N), cmap='viridis', zorder=3, s=15) plt.plot(Y,r,color="black", zorder=3, linewidth=0.9) plt.show() return Y, r par = {'a' : 0.5625, 'Z' : 50, 'Y0': 62, 'psi': 0.4, 'r0': 15, 'N': 150, 'phi':0.2, 'k':0.25, 'b':1.525, 'h':0.5, 'M':8, 'P':1 } Y,r = ISLMb(**par) plt.plot(Y) plt.xlabel('time') plt.ylabel('Y') plt.title('Time series for Income: convergence to the equilibrium') plt.show() plt.plot(r) plt.xlabel('time') plt.ylabel('r') plt.title('Time series for Interest Rate: convergence to the equilibrium') plt.show()

...with Fiscal Policy after T periods

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio ISLM Model case b) with fiscal policy after T-periods """ import numpy as np import matplotlib.pyplot as plt def ISLMb2(N, Y0, r0, a, Z, b, k, h, M, P, phi, psi): #define a vector of zeros Y=np.zeros(N) r=np.zeros(N) #initial condition r[0]=r0 Y[0]=Y0 z=1/(1-a) #draw IS and LM schedule y=np.linspace(0,100,11) r_IS=Z/b-y/(z*b) r_LM=k*y/h-(M/P)/h plt.plot(y,r_IS,color="blue", label="IS0") Z=75 r_IS2=Z/b-y/(z*b) plt.plot(y,r_IS2,color="blue", linestyle='dashed', label="IS1") plt.plot(y,r_LM,color="red", label="LM") plt.legend() plt.grid(True) plt.axhline(color='black', lw=0.5) plt.axvline(color='black', lw=0.5) Z=50 for t in range(55): Y[t+1]=Y[t]+phi*(Z+a*Y[t]-b*r[t]-Y[t]) r[t+1]=r[t]+psi*(k*Y[t]-h*r[t]-(M/P)) Z=75 for t in range(55,N-1): Y[t+1]=Y[t]+phi*(Z+a*Y[t]-b*r[t]-Y[t]) r[t+1]=r[t]+psi*(k*Y[t]-h*r[t]-(M/P)) # plt.plot([0,max(Y)], [0, max(r)], color='black', linestyle='dashed', linewidth=0.5) plt.xlabel('Y') plt.ylabel('r') plt.title('Dynamic ISLM Model B') plt.scatter(Y,r, c=np.arange(N), cmap='viridis', zorder=3, s=15) plt.plot(Y,r,color="black", zorder=3, linewidth=0.9) plt.show() return Y, r par = {'a' : 0.5625, 'Z' : 50, 'Y0': 30, 'psi': 0.4, 'r0': 5, 'N': 150, 'phi':0.2, 'k':0.25, 'b':1.525, 'h':0.5, 'M':8, 'P':1 } Y,r = ISLMb2(**par) plt.plot(Y) plt.xlabel('time') plt.ylabel('Y') plt.title('Time series for Income: convergence to the equilibrium') plt.grid(True) plt.show() plt.plot(r) plt.xlabel('time') plt.ylabel('r') plt.title('Time series for Interest Rate: convergence to the equilibrium') plt.grid(True) plt.show()

4. Lotka Volterra

Jacobian Matrix of Lotka-Volterra (plus Determinant and Trace):

Jacobian Script et al.

Copy to Clipboard
import numpy as np import matplotlib.pyplot as plt import sympy as sym from sympy import * from sympy.matrices import Matrix x=symbols('x') diff=diff(x**2, x) print(diff) def Jacobian(v_str, f_list): vars = sym.symbols(v_str) f = sym.sympify(f_list) J = sym.zeros(len(f),len(vars)) for i, fi in enumerate(f): for j, s in enumerate(vars): J[i,j] = sym.diff(fi, s) return J A=Jacobian('x y', ['(1+r)*x-s*x**2-b*x*y','c*y*x-e*y']) print(A) B=A.det() print(B) C=A.trace() print(C)

TOOLS

Python (project, anaconda)
LateX (project, example)

CompEc 2022/2023

Computational Economics Università degli Studi di Firenze

SYLLABUS

A TENTATIVE PROGRAM

1. Income Expenditure Model

2. Logistic model

3. Dynamic ISLM

4. Lotka Volterra

TOOLS

Computational Economics
Università degli Studi di Firenze