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

Python Scripts:

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

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)

Lotka-Volterra

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[N-1000:N],y[N-1000:N], color='red', s=1) plt.xlabel("prey") plt.ylabel("predator") return x, y parameters = { "x0": 1.5, "y0": 0.5, "r" : 0.5, "s" : 0.21, #0.18,0.16, 0.125 "b" : 1, "e" : 0.6, "c" : 1, "N" : 2000, } x,y= lotka_volterra(**parameters) plt.show() plt.plot(x[0:100], label="prey") plt.plot(y[0:100], label="predator") plt.legend(bbox_to_anchor=(-0.4, 1), loc='upper left', borderaxespad=0.1) 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()

Steady States for the Lotka-Volterra Model

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Thu Mar 16 20:18:54 2023 @author: Giorgio Ricchiuti @website: www.grarchive.net Steady States for the Lotka-Volterra Model """ import sympy as sym #from sympy import * def Lotka_Volterra_SS(): y,x, s,c, r, b, e=sym.symbols('y, x, s, c, r, b, e') system=sym.solve([(1 + r)*x-s*x**2-b*x*y-x, y*c*x-y*e-y], (x,y))#, force=True, manual=True, set=True ) print(system) system=Lotka_Volterra_SS()

5. 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
Python Script:

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 Diagram

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @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)

6. Dynamic effects of increasing heterogeneity in financial markets

Paper:
1. Naimzada A. and Ricchiuti G. (2009) `Dynamic Effects of Increasing Heterogeneity in Financial Markets’, Chaos Solitons and Fractals, Volume 41, Issue 4, August 2009, Pages 1764-1772
2. Naimzada A. and Ricchiuti G. (2012) `A Note on Biased Fundamentalists’, Chaos Solitons and Fractals, Volume 45, Issue 3, March 2012, Pages 224 – 228
Python Script:

Dynamics

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio www.grarchive.net Naimzada A. and Ricchiuti G. (2009) `Dynamic Effects of Increasing Heterogeneity in Financial Markets’, Chaos Solitons and Fractals, Volume 41, Issue 4, August 2009, Pages 1764-1772 https://doi.org/10.1016/j.chaos.2008.07.022 """ import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm def fundamentalist09(N, p0, n0, alfa,beta, gamma, F1,F2): p=np.zeros(N) p[0]=p0 q1=np.zeros(N) q2=np.zeros(N) n1=np.zeros(N) n1[0]=n0 n2=np.zeros(N) n2[0]=1-n1[0] R=np.zeros(N) R[0]=0 z=np.linspace(0,20,501) y=z+beta*alfa*((F2-z)-((F2-F1)/(1+np.exp(gamma*(F2-F1)*(2*z-F1-F2))))) plt.plot(z,y) 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): q1[t]=(F1-p[t]) q2[t]=(F2-p[t]) p[t+1]=p[t]+beta*alfa*((F2-p[t])-((F2-F1)/(1+np.exp(gamma*(F2-F1)*(2*p[t]-F1-F2))))) R[t+1]=(p[t+1]-p[t])/p[t] n1[t+1]=1/(1+(np.exp(gamma*(F1-p[t])**2-gamma*(F2-p[t])**2))) n2[t+1]=1-n1[t+1] plt.plot([p[t],p[t]],[p[t],p[t+1]], color='blue') plt.plot([p[t],p[t+1]],[p[t+1],p[t+1]], color='red') plt.plot([0,max(p)], [0, max(p)], color='black', linestyle='dashed', linewidth=0.5) plt.show() return p , n1, n2, R # par = {'N' : 100, # 'p0' : 6.5, # 'n0': 0.55, # 'alfa': 1.1, # 'beta': 2.1, # 'gamma': 0.3, # 'F1': 8, # 'F2': 15, # } #(N, p0, n0, alfa,beta, gamma, F1,F2) p, n1, n2 , R = fundamentalist09(100,6.5,0.55,1.1,2.1,0.3,8,15) #p, n1, n2 , R = fundamentalist09(100,6.5,0.55,1.1,1,0.8,7,8) plt.plot(p) plt.title("Time series of P") plt.show() plt.plot(n1) plt.plot(n2) plt.title("Time series of Quotas") plt.show() plt.plot(R) plt.title("Time series of Return") plt.show() fig = sm.qqplot(R, line='45') plt.show() a= plt.acorr(R) plt.xlabel("Lags") plt.grid(True) plt.show() plt.hist(R, bins=25) plt.gca().set(title='Frequency Histogram', ylabel='Frequency') plt.show() p, n1, n2 , R = fundamentalist09(100,6.5,0.5,1.1,1,0.8,6,8)

Bifurcation Diagram

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio www.grarchive.net Bifurcation diagram for Naimzada A. and Ricchiuti G. (2009) `Dynamic Effects of Increasing Heterogeneity in Financial Markets’, Chaos Solitons and Fractals, Volume 41, Issue 4, August 2009, Pages 1764-1772 https://doi.org/10.1016/j.chaos.2008.07.022 """ import numpy as np import matplotlib.pyplot as plt def fundamentalist09(N, p0, alfa, beta, gamma, F1,F2): p=np.zeros(N) p[0]=p0 for t in range(N-1): p[t+1]=p[t]+beta*alfa*((F2-p[t])-((F2-F1)/(1+np.exp(gamma*(F2-F1)*(2*p[t]-F1-F2))))) return p # Create the bifurcation diagram def bifurcation_diagram(seed, n_skip, n_iter, step, beta_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=6 F2=8 alfa=1.1 gamma=0.8 p0=seed # Create the A values to loop. For each A value we will plot n_iter points beta_range = np.linspace(beta_min, 5, int(1/step)) for beta in beta_range: x = seed; 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(beta) P.append(x) x = fundamentalist09(N, p0, alfa, beta, gamma, 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(5, 9) plt.xlim(beta_min, 5) plt.xlabel(r'$\beta$') plt.ylabel('P') return A, P, x, X, A2, P2 A, P, p, X, A2, P2 =bifurcation_diagram(6.5, 150, 10, step=0.01, beta_min=3) A, P, p, X, A2, P2 =bifurcation_diagram(7.5, 150, 10, step=0.01, beta_min=3)

7. A Rational Route to Randomness

William A. Brock & Cars H. Hommes, 1997. “A Rational Route to Randomness,” Econometrica, Econometric Society, vol. 65(5), pages 1059-1096, September.
Python Script:

Dynamics

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio www.grarchive.net William A. Brock & Cars H. Hommes, 1997. "A Rational Route to Randomness," Econometrica, Econometric Society, vol. 65(5), pages 1059-1096, September. https://doi.org/10.2307/2171879 """ import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.graphics.tsaplots import plot_acf def BH1997(N, p0, n0, a,beta,s, d,w, C): p=np.zeros(N) #p1e=np.zeros(N) #p2e=np.zeros(N) n1=np.zeros(N) n2=np.zeros(N) #m=np.zeros(N) pi1=np.zeros(N) pi2=np.zeros(N) U1=np.zeros(N) U2=np.zeros(N) R=np.zeros(N) R[0]=0 p_ss=np.zeros(N) p_ss[0]=a/(d+s) p[0]=p0 n1[0]=n0 n2[0]=1-n1[0] pi1[0]=0 pi2[0]=0 U1[0]=0 U2[0]=0 # plt.axhline(color='black', lw=0.5) # plt.axvline(color='black', lw=0.5) for t in range(1, N-1): #update price p[t]=(a-n2[t-1]*s*p[t-1])/(d+n1[t-1]*s) #compute profits pi1[t]=(s/2)*p[t]*p[t]-C pi2[t]=(s/2)*p[t-1]*(2*p[t]-p[t-1]) #compute utility U1[t]=(1-w)*pi1[t]+w*U1[t-1] U2[t]=(1-w)*pi2[t]+w*U2[t-1] #fraction Z= np.exp(beta*U1[t])+np.exp(beta*U2[t]) n1[t]=np.exp(beta*U1[t])/Z n2[t]=1-n1[t] R[t]=(p[t]-p[t-1])/p[t-1] p_ss[t]=p_ss[t-1] return p, pi1, pi2, n1, n2 , U1, U2, R, p_ss par = {'N' : 100000, #time iterations 'p0' : 0.1, 'n0': 0.55, 'a': 3, #demand intercept 'beta': 20, #intensity of choice 'd': 0.5, # demand slope 's': 1.35, #supply slope 'C': 0.2, # cost of the fundamental predictor 'w': 0.1 #fitness measure persistence } p, pi1, pi2, n1, n2 , U1, U2, R, p_ss = BH1997(**par) plt.plot(p) plt.plot(p_ss, color='red') plt.title("Price") plt.xlabel('Time') plt.ylabel('Price') plt.xlim([99000, 99500]) plt.show() plt.plot(n1) plt.plot(n2) plt.xlabel('Time') plt.ylabel('Fraction of rational agents') plt.title("Time series of Quotas") plt.xlim([99900, 99950]) plt.show() plt.plot(U1) plt.plot(U2) plt.xlabel('Time') plt.ylabel('utility') plt.title("Time series of U") plt.xlim([99900, 99950]) plt.show() plot_acf(p) plt.show() fig = sm.qqplot(R, line='45') plt.show() plot_acf(R) plt.show() plt.hist(R, bins=25) plt.gca().set(title='Frequency Histogram', ylabel='Frequency') plt.show() plt.scatter(p,n1, marker='.', alpha=0.1, lw=0.05); plt.xlabel('Price') plt.ylabel('Fraction of rational agents') plt.title('Phase diagrams and homoclinic orbits') plt.show()

Bifurcation Diagram

Copy to Clipboard
#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @author: giorgio Bifurcation Diagram for William A. Brock & Cars H. Hommes, 1997. "A Rational Route to Randomness," Econometrica, Econometric Society, vol. 65(5), pages 1059-1096, September. https://doi.org/10.2307/2171879 """ import numpy as np import matplotlib.pyplot as plt def BH1997(N, p0, n0, a,beta,s, d,w, C): p=np.zeros(N) x=np.zeros(N) n1=np.zeros(N) n2=np.zeros(N) pi1=np.zeros(N) pi2=np.zeros(N) U1=np.zeros(N) U2=np.zeros(N) R=np.zeros(N) R[0]=0 p_ss=np.zeros(N) p_ss[0]=a/(d+s) p[0]=p0 n1[0]=n0 n2[0]=1-n1[0] pi1[0]=0 pi2[0]=0 U1[0]=0 U2[0]=0 for t in range(1, N-1): #update price p[t]=(a-n2[t-1]*s*p[t-1])/(d+n1[t-1]*s) #compute profits pi1[t]=(s/2)*p[t]*p[t]-C pi2[t]=(s/2)*p[t-1]*(2*p[t]-p[t-1]) #compute utility U1[t]=(1-w)*pi1[t]+w*U1[t-1] U2[t]=(1-w)*pi2[t]+w*U2[t-1] #fraction Z= np.exp(beta*U1[t])+np.exp(beta*U2[t]) n1[t]=np.exp(beta*U1[t])/Z n2[t]=1-n1[t] R[t]=(p[t]-p[t-1])/p[t-1] p_ss[t]=p_ss[t-1] x[t]=p[t]-p_ss[t] return p, pi1, pi2, n1, n2 , U1, U2, R, p_ss,x # # Create the bifurcation diagram def bifurcation_diagram(seed, n_skip, n_iter, step, beta_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 N=200 p0=0.1 n0=0.55 a=3 s=1.35 d=0.5 w=0 C=0.5 # Create the A values to loop. For each A value we will plot n_iter points beta_range = np.linspace(beta_min, 6, int(1/step)) for beta in beta_range: p = seed 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): if i >= n_skip: A.append(beta) P.append(p) X.append(x) p, pi1, pi2, n1, n2 , U1, U2, R, p_ss,x=BH1997(N, p0, n0, a,beta,s, d,w, C) # 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) X2=np.delete(X, 1, axis=0) X2=np.delete(X2, np.s_[:100], axis=1) plt.plot(A2, X2, ls='', marker=',', color='red') plt.ylim(-1.5, 1.5) plt.xlim(beta_min, 6) plt.xlabel(r'$\beta$=intensity of choice') plt.ylabel('X') return A, P, X, p, x A, P, X, p, x =bifurcation_diagram(0.5, 150, 100, step=0.01, beta_min=0)

8. Heterogeneous beliefs and routes to chaos in a simple asset pricing model

William A. Brock, Cars H. Hommes, Heterogeneous beliefs and routes to chaos in a simple asset pricing model, Journal of Economic Dynamics and Control, Volume 22, Issues 8–9, 1998, Pages 1235-1274, ISSN 0165-1889,
Python Script:

Dynamics (Two Types)

Copy to Clipboard
""" @author: GR, www.grarchive.net BH1998 with two types of agents """ import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm def BH1998two(N, beta,b1, b2, g1, g2): x=np.zeros(N) # difference between price and optimal price y=np.zeros(N) # delayed difference between price and optimal price n1=np.zeros(N) # quotas of first group n2=np.zeros(N) # quotas of second group z=np.zeros(N) pi1=np.zeros(N) pi2=np.zeros(N) U1=np.zeros(N) U2=np.zeros(N) V=np.zeros(N) s=np.zeros(N) R=1.1 # risk-free C=1 # cost of the fundamental predictor w=0 #fitness measure persistence n10=0.2 delta=0.01 V[0]=0 x[0]=0.21 y[0]=0.2 n1[0]=n10 n2[0]=1-n10 pi1[0]=0 # profits first group pi2[0]=0 # profits second group U1[0]=0 #Utility 1 U2[0]=0 #Utility 2 z[0]=0 # plt.axhline(color='black', lw=0.5) # plt.axvline(color='black', lw=0.5) meanX=np.zeros(N) for t in range(N-1): s[t]=n1[t]*(g1*x[t]+b1)+n2[t]*(g2*x[t]+b2) #update deviation price from fundamental price x[t+1]=s[t]/R #compute profits pi1[t+1]=(x[t+1]-R*x[t])*(g1*y[t]+b1-R*x[t+1])-C pi2[t+1]=(x[t+1]-R*x[t])*(g2*y[t]+b2-R*x[t+1]) #compute utility U1[t+1]=(1-w)*pi1[t+1]+w*U1[t] U2[t+1]=(1-w)*pi2[t+1]+w*U2[t] #fraction z[t+1]= np.exp(beta*U1[t+1])+np.exp(beta*U2[t+1]) n1[t+1]=(delta*n1[t]+(1-delta)*np.exp(beta*U1[t+1]))/z[t+1] n2[t+1]=1-n1[t+1] V[t+1]=(x[t+1]-x[t])/x[t] y[t+1]=x[t] x[t]=x[t+1] return x, y, s, pi1, pi2, n1, n2, U1, U2, V """ Moving Average Function """ #BH1998two(N, beta,b1, b2, g1, g2): x, y, s, pi1, pi2, n1, n2, U1, U2, V = BH1998two(10000, 3.7, b1=0, b2=0, g1=0, g2=1.2) plt.plot(x) plt.title("Deviation from Fundamental Value") plt.xlabel('Time') plt.ylabel('X') plt.xlim([9000, 10000]) plt.ylim([-0.2, 3]) plt.show() plt.plot(n1, label='Fundamentalists') plt.plot(n2, label='Chartists') plt.xlabel('Time') plt.ylabel('Fraction of agents') plt.title("Time series of Quotas") plt.legend() plt.xlim([9450, 9500]) plt.ylim([-0.1, 1.1]) plt.show() fig = sm.qqplot(x, line='45') plt.show() fig = sm.qqplot(V, line='45') plt.show() a= plt.acorr(x) plt.xlabel('Lags') plt.grid(True) plt.show() b= plt.acorr(V) plt.xlabel('Lags') plt.grid(True) plt.show() plt.hist(x, bins=25) plt.xlabel('Deviation from SS') plt.gca().set(title='Frequency Histogram', ylabel='Frequency') plt.show() plt.scatter(x,n1, marker='.', alpha=0.1, lw=0.05); plt.xlabel('Price (deviation from ss)') plt.ylabel('Fraction of rational agents') plt.title('Phase diagrams and homoclinic orbits') plt.show() plt.plot(pi1, label='Profits Fundamentalists') plt.plot(pi2, label='Profits Chartists') plt.xlabel('Time') plt.ylabel('Profits') plt.title("Time series of Profits") plt.legend() plt.xlim([9450, 9500]) plt.ylim([-2,2]) plt.show()

Bifurcation Diagram (Two Types)

9. A Dynamic Exchange Rate Model with Heterogenous Agents

Gori, M., Ricchiuti, G. A dynamic exchange rate model with heterogeneous agents. J Evol Econ 28, 399–415 (2018)
Python Script:

Dynamics

Copy to Clipboard
""" @author: GR, www.grarchive.net Gori and Ricchiuti (2018) `A dynamic exchange rate model with heterogeneous agents' """ import numpy as np import matplotlib.pyplot as plt def gori_ricchiuti(N, E0, Y10, Y20, S0, A1, A2, mpc1, mpc2, mpi1, mpi2, d1, d2, f1, f2,c1, c2, b1, b2, a1, a2, beta, gamma): Y1=np.zeros(N) Y2=np.zeros(N) E=np.zeros(N) S=np.zeros(N) lambda1=np.zeros(N) lambda2=np.zeros(N) R=np.zeros(N) G1=np.zeros(N) G2=np.zeros(N) z1=np.zeros(N) z2=np.zeros(N) k1=np.zeros(N) k2=np.zeros(N) o1=np.zeros(N) o2=np.zeros(N) s=np.zeros(N) g1=np.zeros(N) g2=np.zeros(N) g3=np.zeros(N) Y1[0]=Y10 Y2[0]=Y20 E[0]=E0 S[0]=1/E[0] lambda1[0]=0.4 lambda2[0]=0.3 for t in range(N-1): G1[t]=(A1+mpc1*Y1[t]+(mpi2*Y2[t])/E[t]-mpi1*Y1[t]-Y1[t]) G2[t]=(A2+mpc2*Y2[t]+(mpi1*Y1[t]*E[t])-mpi2*Y2[t]-Y2[t]) z1[t]=(mpi2*Y2[t])/(mpi1*Y1[t]) z2[t]=E[t] S[t]=1/E[t] k1[t]=d1*((2)/(np.exp(-z1[t]+z2[t])+1)-1) k2[t]=f2*((2)/(np.exp(-z1[t]+z2[t])+1)-1) o1[t]=np.exp(-beta*(k1[t])**2) o2[t]=np.exp(-beta*(k2[t])**2) s[t]=(np.exp(-beta*(k1[t])**2)+np.exp(-beta*(k2[t])**2)) lambda1[t+1]=o1[t]/(s[t]) lambda2[t+1]=o2[t]/(s[t]) g1[t]=a2*((2)/(np.exp(-G1[t])+1)-1) g2[t]=b2*((2)/(np.exp(-G2[t])+1)-1) g3[t]=c2*((2)/(np.exp(-lambda1[t+1]*gamma*k1[t]-lambda2[t+1]*gamma*k2[t])+1)-1) Y1[t+1]=Y1[t]+g1[t] Y2[t+1]=Y2[t]+g2[t] E[t+1]=E[t]+g3[t] S[t+1]=1/E[t+1] R[t+1]=(E[t+1]-E[t])/E[t] return Y1, Y2, E, S, R, z1 par = {'N' : 30000, 'E0' : 1.6, 'Y10': 48000, 'Y20': 48000, 'A1': 5000, 'A2': 5000, 'mpc1':0.9, 'mpc2':0.9, 'mpi1':0.1, 'mpi2':0.1, 'a1': 10, 'a2': 10, 'b1': 2, 'b2': 2, 'c1': 2, 'c2': 2, 'd1': 30, 'd2': 2, 'f1': 2, 'f2': 2, 'gamma': 1.05, 'beta': 0.071, 'S0': 1/1.6 } Y1, Y2, E, S, R, z1 =gori_ricchiuti(**par) plt.plot(E) plt.plot(z1) plt.ylim(0,2) plt.xlim([900, 995]) plt.xlabel('time') plt.ylabel('E and Fundamental value') plt.show() plt.hist(E, bins=25, density=True) plt.xlabel('Spot Exchange Rate') plt.gca().set(title='Frequency Histogram', ylabel='Density') plt.show() plt.plot(Y1) plt.plot(Y2) plt.ylim(48500, 51000) plt.xlim([1000, 29950]) plt.xlabel('Time') plt.ylabel('GDP in Y1 and Y2') plt.show()

10. Network: An Introduction (by Francesco Campigli)

Notes by Francesco Campigli
Python script:

An example

Copy to Clipboard
""" Created on March 2023 @author: F. Campigli """ # Import modules import networkx as nx import matplotlib.pyplot as plt ###### Create Graphs G = nx.Graph() # It is an indirect graph ### Other possibilities: DG = nx.DiGraph() #it is an undirected graph # G = nx.MultiGraph() #it is a multi graph # G = nx.MultiDiGraph() #it is a directed multi graph # Indirect graph ###### Create nodes and edges G = nx.Graph([('Marta', 'Luca'),]) G.add_node('Gino') # Add a single node G.add_node('Erica') # Add a single node by mistake G.add_nodes_from(['Levante', 'Libero', 'Elena','Marco','Chiara']) # Add a list of nodes G.add_edge('Marco', 'Marta') # Add one edge, both ends exist G.add_edge('Marta', 'Levante') # Add one edge by mistake G.add_edge('Elena', 'Giulio') # Add one edge, one end does not exist G.add_edges_from([('Levante', 'Elena'), ('Levante', 'Chiara')]) G.add_edge('Gino', 'Massimo') G.add_edge('Chiara', 'Gino') G.add_edge('Gino', 'Libero') G.add_edge('Libero', 'Levante') G.add_edge('Levante', 'Libero') G.add_edge('Libero', 'Luca') G.add_edge('Luca','Chiara') ## Self edge G.add_edge('Chiara','Chiara') ### Example of Direct graph DG.add_edge('Chiara', 'Gino') DG.add_edge('Gino', 'Libero') DG.add_edge('Libero', 'Levante') DG.add_edge('Levante', 'Libero') DG.add_edge('Libero', 'Luca') DG.add_edge('Luca','Chiara') ###### Remove nodes and edges G.remove_node('Giulio') G.remove_node('Erica') G.remove_nodes_from(['Giulio',]) # Safe to remove a missing node using a list G.remove_edge('Luca', 'Chiara') ###### Print information from the graph print(G.nodes()) print(G.edges()) ###### Plot the graph plt.figure(figsize=(6,6)) nx.draw_networkx(G, width = 2) # plt.figure(figsize=(6,6)) # nx.spring_layout(G) # nx.draw_networkx(G) # In[8]: plt.figure(figsize=(6,6)) nx.draw_networkx(DG, width = 2) # In[9]: ###### Create and plot weighted graphs A = nx.Graph() A.add_edge('Giulia', 'Giorgio', weight=25) A.add_edge('Giulia', 'Ugo', weight=2) A.add_edge('Giulia', 'Alessia', weight=2) A.add_edge('Ugo', 'Giorgio', weight=5) A.add_edge('Ugo', 'Alessia', weight=5) A.add_edge('Ugo', 'Lucio', weight=10) A.add_edge('Ugo', 'Francesco', weight=5) A.add_edge('Ugo', 'Francesco', weight=2) plt.figure(figsize=(5,5)) nx.draw_networkx(A) plt.show() # In[10]: plt.figure(figsize=(6,6)) elarge = [(u,v) for (u,v,d) in A.edges(data=True) if d['weight']>=10] emid = [(u,v) for (u,v,d) in A.edges(data=True) if ((d['weight']<10) and (d['weight']>=4))] esmall = [(u,v) for (u,v,d) in A.edges(data=True) if d['weight']<4] pos=nx.circular_layout(A) # positions for all nodes #nodes nx.draw_networkx_nodes(A,pos,node_size=200) #edges nx.draw_networkx_edges(A,pos,edgelist=elarge, width=5) nx.draw_networkx_edges(A,pos,edgelist=emid, width=3,edge_color='b') # , style='dotted' nx.draw_networkx_edges(A,pos,edgelist=esmall, width=2, alpha=0.5, edge_color='r') nx.draw_networkx_labels(A, pos, font_size=20, font_family='sans-serif') plt.axis('off') plt.show() # In[11]: pos = nx.spring_layout(A) betCent = nx.betweenness_centrality(A, normalized=True, endpoints=True) node_color = [20000.0 * A.degree(v) for v in A] node_size = [v * 10000 for v in betCent.values()] plt.figure(figsize=(10,10)) nx.draw_networkx(A, pos=pos, with_labels=True, node_color=node_color, node_size=node_size ) plt.axis('off'); sorted(betCent, key=betCent.get, reverse=True)[:5] plt.show() # In[12]: ### Network metrics print(G.nodes()) # nx.barycenter(G) # Calculate barycenter of a graph. # In[13]: print(G.edges()) # In[14]: # Size L = G.size() print(L) # In[15]: # Number of nodes m = G.number_of_nodes() print(m) # In[16]: # Number of Edges n = G.number_of_edges() print(n) # In[17]: # Density D = 2*m/((n)*(n-1)) print(D) # In[18]: # Degrees nx.degree(G) # In[19]: # Eccentricity nx.eccentricity(G) # In[23]: # Center nx.center(G) # In[24]: # Diameter nx.diameter(G) # In[21]: # Same for Directed graphs L = DG.size() print(L) # In[22]: # Number of Edges n = DG.number_of_edges() print(n) # In[23]: # Number of nodes m = DG.number_of_nodes() print(m) # In[24]: # Density D = m/((n)*(n-1)) print(D) # etc etc # In[25]: nx.degree(DG) # In[ ]: pip

11. Topology Measures

Notes on Topology Measures
Simple Dataset
Python Script:

An example

Copy to Clipboard
""" @author: GR, www.grarchive.net """ import networkx as nx import matplotlib.pyplot as plt import pandas as pd df= pd.read_csv('net_prova1.txt') df.info() df.tail G_fb=nx.read_weighted_edgelist('net_prova1.txt', create_using=nx.DiGraph()) dgr=nx.degree_centrality(G_fb) print(nx.info(G_fb)) print(dgr) plt.figure(figsize=(20,20)) nx.draw_networkx(G_fb)
Python Script 2:

Topology Measurements

Copy to Clipboard
""" @author: www.grarchive.net Network Trade """ import networkx as nx import matplotlib.pyplot as plt import numpy as np import pandas as pd from pathlib import Path import os import glob def trade_network(data_path, year): data = pd.read_csv(data_path) data = data.query("year == @year") G=nx.from_pandas_edgelist(data, 'reportercode', 'partnercode', ['weight2'], create_using=nx.Graph) nx.density(G) print(nx.info(G)) return G src_path = r"FileName.csv" #x = range(2003, 2004) x= [1996,1997] for year in x: G = trade_network('XX.csv', year) dgr=nx.degree_centrality(G) clo=nx.closeness_centrality(G, u=None,distance='weight2') eig=nx.eigenvector_centrality(G, max_iter=100, tol=1e-06, nstart=None, weight='weight2') katz=nx.katz_centrality(G, alpha=0.1, beta=1.0, max_iter=1000, tol=1e-06, nstart=None, normalized=True, weight='weight2') btw=nx.betweenness_centrality(G,k=None, normalized=True, weight='weight2', endpoints=True, seed=None) avd=nx.average_neighbor_degree(G, nodes=None, weight='weight2') ecc=nx.eccentricity(G, v=None, sp=None) clust=nx.clustering(G, weight='weight2') sqclust=nx.square_clustering(G) pagerank=nx.pagerank(G, alpha=0.85, max_iter=100, weight='weight2') centralities=pd.concat([pd.Series(c) for c in (dgr, clo, eig, katz, btw, avd, ecc, clust, sqclust, pagerank)], axis=1 ) centralities.columns=("Degree","Closeness", "Eigenvector", "Katz", "Betweenness", "Av N. Degree", "Eccentricity", "Clustering", "Squared Clustering", "Pagerank") centralities['year']=year os.chdir("new") p = Path(src_path).parent.joinpath(f"{year}_centrality.csv") centralities.to_csv(p) # nx.draw(G, pos=nx.kamada_kawai_layout(G)) # plt.title(year) # plt.show() os.chdir('..') #We will now match the file pattern (‘.csv’) within all of the files located in the current working directory. file_extension = '.csv' all_filenames = [i for i in glob.glob(f"*{file_extension}")] print(f"These are all of the filenames ending in .csv {all_filenames}.") combined_csv_data = pd.concat([pd.read_csv(f, delimiter=',') for f in all_filenames]).to_csv('centrality_complete.csv')

12. Asset price dynamics with heterogeneous beliefs and local network interactions

Mock Exams 1

Text
An Handout of the Solution
Python Script

Copy to Clipboard
""" grarchive.net solution Mock 1 """ import numpy as np import matplotlib.pyplot as plt import sympy as sym from sympy import * from sympy.matrices import Matrix def funMock1(): y,x,g,b =sym.symbols('y, x, g, b') system=sym.solve([1-g*(x**2)+y-x, b*x-y], (x,y))#, force=True, manual=True, set=True ) print(system) system=funMock1() 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-g*(x**2)+y','b*x']) print(A) B=A.det() print(B) C=A.trace() print(C) def Mock1(N, y0, x0, g,b): x=np.zeros(N) y=np.zeros(N) y[0]=y0 #initial condition x[0]=x0 #initial condition for t in range(N-1): x[t+1]=1-g*x[t]**2+y[t] y[t+1]=b*x[t] return x, y x,y= Mock1(N=100, y0=0.5, x0=0.6, g=0.2, b=0.3) N=100 plt.scatter(x[N-1000:N],y[N-1000:N], color='red', s=1) plt.xlabel(r"$x_t$") plt.ylabel(r"$y_t$") plt.show() plt.show() plt.plot(x[0:100], label="$x_t$") plt.plot(y[0:100], label="$y_t$") plt.legend(bbox_to_anchor=(-0.4, 1), loc='upper left', borderaxespad=0.1) plt.show() plt.hist(x, 50, density=True, facecolor='b', alpha=0.75) plt.ylabel('Probability') plt.title('Histogram of $x_t$') plt.grid(True) plt.show() plt.hist(y, 50, density=True, facecolor='g', alpha=0.75) plt.ylabel('Probability') plt.title('Histogram of $y_t$') plt.grid(True) plt.show() def bifurcation(seed, n_skip, n_iter, step, b_min, b_max): y0=0.5 x0=0.6 g=0.2 print("Starting with x0 seed {0}, skip plotting first {1} iterations, then plot next {2} iterations.".format(seed, n_skip, n_iter)); # Array of b values, the x axis of the bifurcation plot B = [] # Array of X_t values, the y axis of the bifurcation plot X1 = [] N=n_skip+n_iter b_range = np.linspace(b_min, b_max, int(1/step)) for b in b_range: x = seed; for i in range(n_iter+n_skip+1): x,y= Mock1(N, y0, x0, g, b) if i >= n_skip: B.append(b) X1.append(x) X12=np.delete(X1, 1, axis=0) X12=np.delete(X12, np.s_[:100], axis=1) B2=np.delete(B, 1, axis=0) plt.plot(B2, X12, ls='', marker=',', color='red') plt.ylim(-3, 3) plt.xlim(b_min, b_max) plt.xlabel(r'$\beta$') plt.ylabel(r'$x_t$') plt.show() return B, X1, X12, B2 B, X1, X12, B2= bifurcation(seed=0.25, n_skip=80, n_iter=120, step=0.01, b_min=0.1, b_max=0.99)

Mock Exams 2

Text
An Handout of the Solution
Python Script

Copy to Clipboard
""" @author: Giorgio Ricchiuti @website: www.grarchive.net Description: solution Mock 2 """ import matplotlib.pyplot as plt import networkx as nx from pathlib import Path import pandas as pd A=nx.DiGraph() A.add_edge('Lapo', 'Nicola', weight=3) A.add_edge('Lapo', 'Giorgio', weight=2) A.add_edge('Giorgio', 'Francesca', weight=3) A.add_edge('Yuri', 'Francesca', weight=4) A.add_edge('Lapo', 'Francesca', weight=4) A.add_edge('Yuri', 'Lapo', weight=1) A.add_edge('Yuri', 'Nicola', weight=1) A.add_edge('Yuri', 'Emma', weight=2) A.add_edge('Emma', 'Nicola', weight=3) A.add_edge('Emma', 'Alice', weight=4) A.add_edge('Francesca', 'Alice', weight=5) A.add_edge('Alice', 'Emma', weight=2) A.add_edge('Alice', 'Giorgio', weight=1) plt.figure(figsize=(6,6)) nx.draw_networkx(A) plt.show() pos = nx.spring_layout(A) betCent = nx.betweenness_centrality(A, normalized=True, endpoints=True) node_color = [20000.0 * A.degree(v) for v in A] node_size = [v * 10000 for v in betCent.values()] plt.figure(figsize=(10,10)) nx.draw_networkx(A, pos=pos, with_labels=True, node_color=node_color, node_size=node_size ) plt.axis('off'); dgr=nx.degree_centrality(A) clo=nx.closeness_centrality(A, u=None, distance='weight') eig=nx.eigenvector_centrality(A, max_iter=100, tol=1e-06, nstart=None, weight='weight') btw=nx.betweenness_centrality(A, k=None, normalized=True, weight='weight', endpoints=True, seed=None) avd=nx.average_neighbor_degree(A, nodes=None, weight='weight') clust=nx.clustering(A, weight='weight') sqclust=nx.square_clustering(A) centralities=pd.concat([pd.Series(c) for c in (dgr, clo, eig, btw, avd, clust, sqclust)], axis=1) centralities.columns=("Degree","Closeness", "Eigenvector", "Betweenness", "Av N. Degree", "Clustering", "Squared Clustering") src_path=r"Filename.csv" p=Path(src_path).parent.joinpath(f"centrality.csv") centralities.to_csv(p)

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

5. Heterogeneous Fundamentalists and Imitative Processes

6. Dynamic effects of increasing heterogeneity in financial markets

7. A Rational Route to Randomness

8. Heterogeneous beliefs and routes to chaos in a simple asset pricing model

9. A Dynamic Exchange Rate Model with Heterogenous Agents

10. Network: An Introduction (by Francesco Campigli)

11. Topology Measures

12. Asset price dynamics with heterogeneous beliefs and local network interactions

Mock Exams 1

Mock Exams 2

TOOLS

Computational Economics
Università degli Studi di Firenze