LabEC 2025/26

Introduzione all’Economia Computazionale

Funzione Produzione Cobb-Douglas

Script

Copy to Clipboard

Note sul Modello Reddito-Spesa (caso statico e caso dinamico)

"""
Dynamic Income–Expenditure model (Keynesian cross).
Parameters:
    N       = number of periods
    x0      = initial income
    A       = autonomous expenditure
    epsilon = c(1-t)-q
    phi     = speed of adjustment
Returns:
    Array with the simulated income trajectory.
"""

import numpy as np
import matplotlib.pyplot as plt

def IE(N, x0, A, epsilon, phi):

# Initialize vector of incomes with zeros
    x = np.zeros(N)
    # Set initial condition
    x[0] = x0
    
    # Create grid for the static IE diagram
    Y = np.linspace(0, 50000, 11)
    Z = A + epsilon * Y
    
    # Plot expenditure line Z(Y) = A + epsilon*Y
    plt.plot(Y, Z, color="black", label=r"$Z(Y)$")
    # Plot the 45-degree line
    plt.plot(Y, Y, label="45° line")
    # Axes crossing at zero
    plt.axhline(color='black', lw=0.5)
    plt.axvline(color='black', lw=0.5)
    
    # Cobweb staircase: dynamic adjustment
    for t in range(N-1):
        # Update rule: Y_{t+1} = Y_t + phi*(Z(Y_t) - Y_t)
        x[t+1] = x[t] + phi*(A + epsilon*x[t] - x[t])
        
        # Vertical segment: (Y_t, Y_t) -> (Y_t, Y_{t+1})
        plt.plot([x[t], x[t]], [x[t], x[t+1]], color='blue')
        # Horizontal segment: (Y_t, Y_{t+1}) -> (Y_{t+1}, Y_{t+1})
        plt.plot([x[t], x[t+1]], [x[t+1], x[t+1]], color='red')
    
    plt.title("Dynamic Income–Expenditure Model")
    plt.xlabel(r"$Y_t$")
    plt.ylabel(r"$Z(Y_t)$")
    plt.legend()
    plt.grid(True)
    plt.show()
    
    # Time series of income
    plt.plot(x, marker="o")
    plt.title("Income time series")
    plt.xlabel("t")
    plt.ylabel(r"$Y_t$")
    plt.grid(True)
    plt.show()
    
    return x

# Example simulation
x = IE(20, 10000, 15000, 0.53, 2.5)

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

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)

Modello Gradiente (articolo)

Dynamics

Copy to Clipboard

Bifurcation

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)

Modello Agenti Eterogenei con Due Fondamentalisti (articolo)

Dynamics

Copy to Clipboard

# -*- coding: utf-8 -*-
"""
Created on Fri Mar 10 13:19:01 2023

@author: Giorgio
Heterogeneous Fundamentalists and Imitative Processes
Naimzada A. and Ricchiuti G. (2008)
Applied Mathematics and Computation, Volume 199, Issue 1, 15 May 2008, Pages 171-180.
https://doi.org/10.1016/j.amc.2007.09.061
 
"""

import numpy as np 
import matplotlib.pyplot as plt

def fundamentalist(N, p0, a, F1, F2):
    p=np.zeros(N)
    p[0]=p0
    n1=np.zeros(N)
    n2=np.zeros(N)
    n1[0]=0.4
    n2[0]=1-n1[0]
    
    z=np.linspace(0,5,501)
    y=z+a*((F1-z)*(F2-z)*(F1+F2-2*z)/((F1-z)**2+(F2-z)**2))
    plt.plot(z,y,color='black')
    plt.plot(z,z)
    plt.grid(True)
    plt.axhline(color='black', lw=0.5)
    plt.axvline(color='black', lw=0.5)
    
    for t in range(N-1):
        p[t+1]=p[t]+a*((F1-p[t])*(F2-p[t])*(F1+F2-2*p[t])/((F1-p[t])**2+(F2-p[t])**2))
        plt.plot([p[t],p[t]],[p[t],p[t+1]])
        plt.plot([p[t],p[t+1]],[p[t+1],p[t+1]])
        n1[t+1]=1-((F1-p[t])**2)/((F1-p[t])**2+(F2-p[t])**2)
        n2[t+1]=1-n1[t+1]
    return p, n1,n2

"""
Figures
"""

fig=plt.figure()
ax1=fig.add_subplot(321)
p, n1, n2=fundamentalist(100, 3.6, 1, 1, 4)
ax2=fig.add_subplot(322)
plt.plot(p)

ax1=fig.add_subplot(323)
p, n1, n2=fundamentalist(100, 3.6, 2.1, 1, 4)
ax2=fig.add_subplot(324)
plt.plot(p)

ax1=fig.add_subplot(325)
p, n1, n2=fundamentalist(100, 3.6, 2.8, 1, 4)
ax2=fig.add_subplot(326)
plt.plot(p)
plt.show()

p, n1, n2=fundamentalist(100, 3.6, 3.1, 1, 4)
plt.show()
plt.plot(p)
plt.show()

p, n1, n2=fundamentalist(100, 3.6, 3.1, 1, 4)
plt.show()
fig=plt.figure()
ax1=fig.add_subplot(321)
ax1=fig.add_subplot(211)
plt.plot(n1)
plt.ylabel("quotas")
plt.plot(n2)
ax1=fig.add_subplot(212)
plt.plot(p)
plt.ylabel("price")
plt.xlabel("time")

plt.show()

Bifurcation

Copy to Clipboard

"""
@author: giorgio

"""

import numpy as np
import matplotlib.pyplot as plt

def fundamentalist(N,  p0, a,F1,F2):
    p=np.zeros(N)
    p[0]=p0
    for t in range(N-1):
        p[t+1]=p[t]+a*((F1-p[t])*(F2-p[t])*(F1+F2-2*p[t]))/((F1-p[t])**2+(F2-p[t])**2)
        
    return p
    
    #    return p+a*((F1-p)*(F2-p)*(F1+F2-2*p))/((F1-p)**2+(F2-p)**2)

# Create the bifurcation diagram

def bifurcation_diagram(seed, n_skip, n_iter, step, a_min):
    print("Starting with p0 seed {0}, skip plotting first {1} iterations, then plot next {2} iterations.".format(seed, n_skip, n_iter));
    # Array of A values, the x axis of the bifurcation plot
    A = []
    # Array of P_t values, the y axis of the bifurcation plot
    P = []
    X=[]
    N=n_skip+n_iter
    F1=1
    F2=5
    p0=0.5
    # Create the A values to loop. For each A value we will plot n_iter points
    a_range = np.linspace(a_min, 3, int(1/step))

for a in a_range:
        x = seed;
   
        # For each A, iterate the function and collect datapoint if n_skip iterations have occurred
        for i in range(n_iter+n_skip+1):
#the skip in this version does not work, I spent all day...then I bypassed this issue
# with lines 52-54. 
            if i >= n_skip:
                A.append(a)
                P.append(x)
              
            x = fundamentalist(N, p0, a,F1,F2)
    # the next three lines, are able to clean 
    A2=np.delete(A, 1, axis=0) 
    P2=np.delete(P, 1, axis=0)
    P2=np.delete(P2, np.s_[:100], axis=1)
    plt.plot(A2, P2, ls='', marker=',', color='red')
    plt.ylim(0, 7)
    plt.xlim(a_min, 3)
    plt.xlabel('alfa*beta')
    plt.ylabel('P')

return A, P, x, X, A2, P2

A, P, p, X, A2, P2 =bifurcation_diagram(0.5, 150, 10, step=0.001, a_min=1.5)

Modello ISLM Dinamico (note)

Casa 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()

Caso 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()

Dynamics 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-periods

Copy to Clipboard

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

LabEC 2025/26

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

PROGRAMMA

SCRIPT E NOTE

Script

Script

Script

Bifurcation py 1

Dynamics

Bifurcation

Dynamics

Bifurcation

Casa a)

Caso b)

Dynamics with fiscal policy

Dynamics with fiscal policy after T-periods

Script

SOFTWARE

LabEC 2025/26

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

PROGRAMMA

SCRIPT E NOTE

SOFTWARE

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