LabEC 2022/23

Cobb Douglas (notes)

Script

Copy to Clipboard

Income Expenditure Model (notes)

Script

Copy to Clipboard

"""
Created on Wed Oct 12 16:19:07 2022
@author: Giorgio Ricchiuti 
@website: www.grarchive.net

Description
Income-Expenditure Model (static and dynamical)
see notes https://grarchive.net/wp-content/uploads/2022/10/NotesIE.pdf

"""

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

def static_IE(Z,z):
    y=sym.Symbol('y')
    eqn=sym.solveset(Z+z*y-y,y)
    print(eqn)

Z=1000
z=0.8

static_IE(Z, z)
controllo=Z/(1-z)
print('Controllo: {}'.format(controllo))
  
def IE(N, x0,Z,z,phi):
    
    x=np.zeros(N)
    x[0]=x0
    y=np.linspace(0, 50000, 11)
    R=Z+z*y
    plt.plot(y,R, color="black")
    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*(Z+z*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.title("Dynamic IE Model")
    plt.show()
    return x

x = IE(20,10000,15000,0.53,1)
plt.plot(x)
plt.title("Time series for Incomce")
plt.show

Logistic Function (articolo Robert May, Nature 1976)

script

Copy to Clipboard

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, 1.3)
plt.title('Logistic Map 2')
plt.grid(True)
plt.show()

h= logistic_map(35, 0.3, 3.2)
plt.title('Logistic Map 3')
plt.grid(True)
plt.show()

plt.plot(f, color='blue')
plt.plot(g, color='red')
plt.plot(h, color='green')
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()

Modello Gradiente (articolo)

Dynamics

Copy to Clipboard

Bifurcation Diagram

Copy to Clipboard

import numpy as np
import matplotlib.pyplot as plt

def cubica(x,a,b,c,k):
    return x+k*(a-c-4*b*x**3)

# Create the bifurcation diagram
def bifurcation_diagram(seed, n_skip, n_iter, step, k_min, a,b,c):
    print("Starting with x0 seed {0}, skip plotting first {1} iterations, then plot next {2} iterations.".format(seed, n_skip, n_iter));
    # Array of K values, the x axis of the bifurcation plot
    K = []
    # Array of x_t values, the y axis of the bifurcation plot
    X = []
    
    # Create the K values to loop. For each K value we will plot n_iter points
    k_range = np.linspace(k_min, 0.29, int(1/step))

for k in k_range:
        x = seed
        # For each K, iterate the cubica and collect datapoint if n_skip iterations have occurred
        for i in range(n_iter+n_skip+1):
            if i >= n_skip:
                K.append(k)
                X.append(x)
                
            x = cubica(x,a,b,c,k);
    # Plot the data    
    plt.plot(K, X, ls='', marker=',')
    plt.ylim(0, 2)
    plt.xlim(k_min, 0.29)
    plt.xlabel('k')
    plt.ylabel('X')
    plt.show()

bifurcation_diagram(1.1, 100, 5, 0.0001, 0.15, 4, 0.6, 0.5)

Dynamic ISLM (note)

Dynamics a)

Copy to Clipboard

"""
@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()

Dynamics 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.8,
      'r0': 12,
      'N': 150,
      'phi':0.05,
      '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()

Dynamic with fiscal policy

Copy to Clipboard

"""
@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.8,
      'r0': 15,
      'N': 150,
      'phi':0.05,
      '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()

Dynamics with fiscal policy after T-period

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.8,
      'r0': 5,
      'N': 150,
      'phi':0.05,
      '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()

Lotka Volterra Model

Script

Copy to Clipboard

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 26 22:03:13 2020

@author: gr

Model Lotka-Volterra 
X_{t+1}= (1+r) X_t-sX^2_{t}-bX_{t}Y_{t}\\
Y_{t+1}= Y_{t}(cX_t-e)

"""

import numpy as np 
import matplotlib.pyplot as plt

def lotka_volterra(N, y0, x0, r, s, b, e, c):

x=np.zeros(N) #prey
    y=np.zeros(N) #predator

y[0]=y0 #initial condition
    x[0]=x0 #initial condition

for t in range(N-1):
        x[t+1]=(1+r)*x[t]-s*(x[t]**2)-b*x[t]*y[t]
        y[t+1]=y[t]*(c*x[t]-e)
        
        
    plt.scatter(x,y, color='red', s=1)    
    plt.xlabel("prey")
    plt.ylabel("predator")
    plt.legend()
    return x, y

parameters = {
    "x0": 1,
    "y0": 0.5,
    "r" : 0.35,
    "s" : 0.125,
    "b" : 1,
    "e" : 0.6,
    "c" : 1,
    "N" : 200,
}
x,y= lotka_volterra(**parameters)
plt.show()
plt.plot(x, label="prey")
plt.plot(y, label="predator") 
plt.legend()
plt.show()
plt.hist(x,  50, density=True, facecolor='b', alpha=0.75)
plt.ylabel('Probability')
plt.title('Histogram of Prey')
plt.grid(True)
plt.show()
plt.hist(y, 50, density=True, facecolor='g', alpha=0.75)
plt.ylabel('Probability')
plt.title('Histogram of Predator')
plt.grid(True)
plt.show()

SIRD Model

Script Discrete SIR(D) plus reinfections

Copy to Clipboard

# -*- coding: utf-8 -*-
"""
Created on Wed Nov 23 17:23:13 2022

@author: GR

SIRD Model 
"""

import numpy as np 
import matplotlib.pyplot as plt

def SIRD(n,S0,I0,sigma, gamma, mu):
    
    S=np.zeros(n) # susceptibles
    I=np.zeros(n) # infected
    R=np.zeros(n) # recovered
    D=np.zeros(n) # deaths
    N=np.zeros(n)  #population
    cumD=np.zeros(n)
    S[0]=S0
    I[0]=I0
    R[0]=0
    D[0]=mu*I0
    N[0]=S[0]+I[0]+R[0]
    cumD[0]=0
    for t in range(n-1):
        
        S[t+1]=S[t]*(1-sigma*I[t]/N[t])
        I[t+1]=I[t]*(1+sigma*S[t]/N[t]-gamma)-D[t]
        R[t+1]=R[t]+gamma*I[t]-mu*I[t]
        N[t+1]=S[t+1]+I[t+1]+R[t+1]
        D[t+1]=mu*I[t]
        cumD[t+1]=cumD[t]+D[t+1]
    
    plt.plot(S, color='red', label="suscetibles")
    plt.plot(I, color='blue', label="infected")
    plt.plot(R, color='green', label="recovered")
    plt.plot(cumD, color='black', label="cum deaths")
    plt.plot(N, color='orange', label="population")
    plt.legend(bbox_to_anchor=(1.05, 1),
                         loc='upper left', borderaxespad=0)
    plt.grid(True)
    plt.xlabel("days")
    plt.ylabel("number of people")    
    
    return S, I, R, D, N, cumD
    
    
S, I, R, D, N, cumD=SIRD(100,1000,2,1.1, 0.18, 0.009)

LabEC 2022/23

Laboratorio di Economia Computazionale
Università degli Studi di Firenze
Dipartimento di Scienze per l’Economia e l’Impresa

PROGRAMMA

SOFTWARE

Script

Script

script

Dynamics

Bifurcation Diagram

Dynamics a)

Dynamics b)

Dynamic with fiscal policy

Dynamics with fiscal policy after T-period

Script

Script Discrete SIR(D) plus reinfections

LabEC 2022/23

Laboratorio di Economia Computazionale Università degli Studi di Firenze Dipartimento di Scienze per l’Economia e l’Impresa

PROGRAMMA

SOFTWARE

Laboratorio di Economia Computazionale
Università degli Studi di Firenze
Dipartimento di Scienze per l’Economia e l’Impresa