In the following, you will find the Python Code to my book on Public Economics (2019). You can also download the code from the download page
https://www.uni-augsburg.de/de/fakultaet/wiwi/prof/vwl/heer/pubec-buch/computer-code/
If you have not installed PYTHON yet or have little prior experience in PYTHON programming, an ideal starting point is provided by the web page of THOMAS J. SARGENT and JOHN STACHURSKI on 'Python Programming for Economics and Finance':
https://python-programming.quantecon.org/intro.html
The order of the following computer codes is chosen with intent. The programs get more sophisticated in each step. We start with applications from the Chapter 3 on overlapping generations model and continue to cover material from Chapter 6 on pensions. We cover the following numerical tools in these exercises:
- functions
- loops
- non-linear equations solver in one variable
- searching for a good initial value as input into the non-linear eqs solver
- non-linear equations solver in multiple variables
- solution of a two-point boundary problem for a first-order non-linear difference equation: Reverse Shooting
- nested application of the non-linear eqs solver fsolve
The numerical and economic problems are described at length in Chapters 3 and 6 of the book Public Economics. Slides, youtube tutorials and solutions to the problem section are also available from the web page of the book (together with MATLAB and Gauss code):
https://www.uni-augsburg.de/de/fakultaet/wiwi/prof/vwl/heer/pubec-buch/
If you have questions or comments, please do not hesitate to contact me: Burkhard.Heer@wiwi.uni-augsburg.de
The computer program Ch3_olg_dyn.py is a good starting point if you have little programming experience. It computes the transition in a simple 2-period OLG model as described in Section 3.2 of my book. The economic problem is described in more detail in the first 5 minutes of my youtube video:
https://www.youtube.com/watch?v=hZCpyQFk5o4&list=PL-anmMgSYtuEHGjWmUUPl8PCwalcjHvpn&index=3
In the application Ch3_olg_dyn.py, we learn the use of 'functions' in Python.
In order to run the program, you need to install Python and the libraries numpy and matlibplot (as described on the pages by Sargent and Stachurski). You can also download the code of the Python code from my webpage:
https://www.uni-augsburg.de/de/fakultaet/wiwi/prof/vwl/heer/pubec-buch/computer-code/
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (10,6)
# parameter values
alpha = 0.36
beta = 0.40
n = 0.1
t0 = 20 # number of transition periods
# steady-state capital stock
k = beta / (1+beta) * ((1-alpha) / (1+n) )**(1/(1-alpha))
# initial capital stock
k0 = k / 3
def kdyn(x):
y = beta / (1+beta) * (1-alpha) / (1+n) * x**(alpha)
return y
kt = [] # initial value of vector
kt.append(k0) # kt[0] = k0
# iterate over periods i=0,..t0-1
for i in range(int(t0)):
k1 = kdyn(k0)
kt.append(k1)
k0=k1
print("rows in kt: " +str(len(kt)))
plt.plot(kt)
plt.show()
rows in kt: 21
The figure displays the transition dynamics in the OLG model and is equivalent to Fig. 3.6 on page 76 in my book on 'Public Economics'.
In the next application, we look at the so-called turnpike behavior as described in Chapter 3.3.4 on 'Dynamics in the Command Optimum' in my book on 'Public Economics'. The economic problem is also described in the following youtube video:
https://www.youtube.com/watch?v=a73hFit4OUQ&list=PL-anmMgSYtuEHGjWmUUPl8PCwalcjHvpn&index=5
The basic new element in this code is the non-linear equations solver 'broyden1' which is provided by the library 'optimize'. The code of Ch3_turnpike.py which you can also download from my homepage is the following:
import scipy.optimize
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (10,6)
# parameter values
alpha = 0.36
beta = 0.40
n = 0.1
bigt = 20 # number of transition periods
k = (alpha/n)**(1/(1-alpha)) # steady state value of capital stock
k0 = 0.7 * k # initial capital stock
kfinal = 0.7 * k # final capital stock
kt = np.zeros(bigt+1) # time path for capital stock
kt[0] = k0
# function that computes the transition given
# initial capital stock k0 and guess k1
# input: x, value of capital stock in period 1
# output: k[bigt]-kfinal
def kdyn(x):
kt[1] = x
for i in range(1,bigt):
temp = (1+alpha*kt[i]**(alpha-1)) / (1+n)**2 * (kt[i-1]+kt[i-1]**alpha - (1+n)*kt[i])
kt[i+1] = (kt[i]+kt[i]**alpha) / (1+n) -temp;
# print(i, kt[i])
return kt[bigt] - kfinal
# tests if kdyn is computing the transition for our guess
kguess = 5.35
ksolution = kdyn(kguess)
x = scipy.optimize.broyden1(kdyn, kguess, f_tol=1e-10)
print("solution for kt[1]: " + str(x))
plt.plot(kt)
plt.show()
solution for kt[1]: 5.336927684600741
In the first Python application of this chapter, we compute the welfare effects and the transition dynamics following the abolition of a pay-as-you-go pension system. The dynamic problem is described in more detail in Chapter 6.3 of my book on 'Public Economics'. In addition, you can find the description of the economic and computational problem in the following youtube video:
https://www.youtube.com/watch?v=SGIxL6IRtHE&list=PL-anmMgSYtuEHGjWmUUPl8PCwalcjHvpn&index=9
The Python code 'Ch6_social_security1.py' can also be downloaded from my homepage:
https://www.uni-augsburg.de/de/fakultaet/wiwi/prof/vwl/heer/pubec-buch/computer-code/
In this program, we, again, have to solve a non-linear equation problem. Therefore, we will make use of the command 'fsolve' from the library optimise. The command 'fsolve' is very convenient. You can use it for one- and multi-dimensional non-linear equations problems.
As an additional new element, we have to search for a good initial value to this non-linear equation problem. First, we have to reformulate the non-linear equation so that the non-linear equations solver does not get stuck at the trivial steady state with zero capital stock (therefore, we divide the non-linear equation by the capital stock). Second, we have to search for a good initial value. For this reason, we specifiy an equispaced grid over the capital stock, compute the respective values of the non-linear equation and allocate the minimum of its absolute values.
Note: The code is not written to be elegant and slim, but rather intuitive for the beginner. Therefore, I did not avoid repitition of code or so-called 'magic numbers'. To improve the code and apply best practises, please consult Chapter 15 of Sargent's and Starchurski's tutorial referenced above.
import scipy.optimize
import numpy as np
import matplotlib.pyplot as plt
import math
plt.rcParams['figure.figsize'] = (10,6)
log = math.log
e = math.e
# part 2: parameterization
labor = 0.30 # constant labor supply
alpha = 0.36 # production elasticity of capital
b = 0.40 # discount factor
n = 0.10 # population growth rate
nu0 = 257.15 # disutility of labor
nu1 = 3.33 # inverse of Frisch elasticity
tau = (0, 0.30) # parameter values for pension replacement rate (tuple)
nt=10 # number of transition periods
# part 3: defining of the functions
#
# wage function
def wage(k,l):
return (1-alpha) * k**alpha * l**(-alpha)
# interest rate function
def interest_rate(k,l):
return alpha * k**(alpha - 1) * l**(1-alpha)
# steady state computation
# input: guess for k
# output: steady-state condition (=0 in steady state)
def ksteady(k):
w = wage(k,labor) # wage
r = interest_rate(k,labor) # interest rate
d = tau0 * w * labor # pension
# non-linear equation (6.19) in Heer (2019)
y = (1+n)*k-w*labor*b/(1+b) + 1/(1+b)*(1+b+b*r+n)/(1+r)*d
# reformulatoin of (6.19) so that k=0 is not a solution
# we are looking for non-trivial steady state k>0
y= y/k
return y
# function that computes k_t+1 given k_t
def kdyn(k1):
w0 = wage(k0,labor) # wage in t
w1 = wage(k1,labor) # wage in t+1
r1 = interest_rate(k1,labor) # interest rate in t+1
d1 = tau1 * w1 * labor # pension in t+1
y = (1+n)*k1 - (1-tau0)*w0*labor*b/(1+b) + 1/(1+b)*(1+n)/(1+r1) * d1
return y
# part 4: main program
tau0 = tau[0] # initial steady state: no pension
# 4.1
# come up with an initial guess for steady state value
# by searching over a grid of values for k
kgrid = np.linspace(0.001, 1, 100)
ygrid = np.zeros(100) # values of non-linear eq for k
for i in range(100):
ygrid[i] = ksteady(kgrid[i])
plt.plot(kgrid, ygrid, 'b-', linewidth=2)
plt.show()
ygrid = abs(ygrid)
# print("Minimum of ygrid: " + str(ygrid.min()))
i0 = np.where(ygrid == ygrid.min()) # find index of y with minimum
kguess = kgrid[i0]
# print("Index of Minimum: " + str(kguess))
# print("ksteady for initial guess: " + str(ksteady(kguess)))
# 4.2
# compute initial steady state with tau=0 and steady-state lifetime utility
kss = scipy.optimize.fsolve(ksteady, kguess)
# print(str(kss))
wss = wage(kss,labor)
rss = interest_rate(kss, labor)
dss = tau0 * wss * labor
c1ss = 1/(1+b) * (wss*labor + (n-rss)/(1+rss)*dss)
c2ss = b*c1ss*(1+rss)
utilss = log(c1ss) + b*log(c2ss)-nu0*labor**(1+nu1)/(1+nu1)
# print(utilss)
# 4.3
# compute final steady state with tau=0.3 and steady-state lifetime utility
tau0 = tau[1]
# find new initial value
for i in range(100):
ygrid[i] = ksteady(kgrid[i])
ygrid = abs(ygrid)
i0 = np.where(ygrid == ygrid.min()) # find index of y with minimum
kguess = kgrid[i0]
kssd = scipy.optimize.fsolve(ksteady, kguess)
# print(ksteady(kssd))
# print(str(kssd))
wssd = wage(kssd,labor)
rssd = interest_rate(kssd, labor)
dssd = tau0 * wssd * labor
c1ssd = 1/(1+b) * (wssd*labor + (n-rssd)/(1+rssd)*dssd)
c2ssd = b*c1ssd*(1+rssd)
utilssd = log(c1ssd) + b*log(c2ssd)-nu0*labor**(1+nu1)/(1+nu1)
# print(utilssd)
# 4.4
# compute steady state welfare effect
cec = e**((utilssd-utilss)/(1+b) )-1
print("Welfare effect in steady state from the introduction of a pension:")
print("consumption equivalent change: " + str(cec))
# 4.5
# computation of the dynamics
# from initial to final steady state
kt = np.zeros(nt+1) # time path for capital stok in periot t
utilt = np.zeros(nt) # lifetime utility of generation born in t
periods = range(nt) # range of period 0,..., nt
cect = np.zeros(nt) # welfare effects, % of consumption
kt[0] = kssd # initial value of capital stock in period 0
utilt[0] = utilssd
cect[0] = e**((utilt[1]-utilss)/(1+b) )-1
tau0 = 0.3 # tax rate at young age
tau1 = 0.3 # tax rate at old age
# i=0: Period 0 with steady state, tau=30% for young and old
# i=1: Period 1, young generation still has to finance old agents
# (tau0=30% so that d_1=tau_0 w_1 l), but will not receive
# a pension in old age in period 2 (i=2) with tau1=0% so that d_2=0
for i in periods:
k0=kt[i]
# print(i,k0)
if i==2:
tau1 = 0 # first generation still has to pay contributions
if i==3:
tau0 = 0 # first generation which does not have to pay contributions
k1 = scipy.optimize.fsolve(kdyn, k0)
kt[i+1] = k1
w0 = wage(k0,labor)
w1 = wage(k1,labor)
r1 = interest_rate(k1, labor)
d0 = tau0 * w0 * labor
d1 = tau1 * w1 * labor
c1 = 1/(1+b)*(w0*labor-d0+d1*(1+n)/(1+r1))
c2 = b * c1 * (1+r1)
utilt[i] = log(c1) + b*log(c2) - nu0*labor**(1+nu1) / (1+nu1)
cect[i] = e**((utilt[i]-utilss) / (1+b) ) -1
fig, axes = plt.subplots(2, 1, figsize=(8, 16))
axes[0].set_xlabel('time')
axes[0].set_ylabel('capital')
axes[0].plot(kt[1:nt])
axes[1].set_xlabel('time')
axes[1].set_ylabel('welfare')
axes[1].plot(cect[1:nt])
plt.show()
Welfare effect in steady state from the introduction of a pension: consumption equivalent change: -0.38428069240593954
The output of the program is as follows: In the top figure, the non-linear steady-state equation (6.19) is pictured over a range of the steady-state capital stock. The absolute minimum of this function provides us with a guess for the solution of the non-linear equation.
Next, the steady state welfare effect of an introduction of a public pay-as-you-go pension system is computed. In steady state, welfare falls by 38% of total consumption.
In the middle figure, the transition of the capital stock is illustrated for the case that pensions are abolished in period 1.
In the bottom figure, the generational welfare during the transition is illustrated. Notice that welfare drops for the generation born in period 1. They still have to finance the pensions of the old generation but do not receive a pension themselves.
In the following example, we look at a 2-period OLG model with elastic labor and study the effects of the abolition of public PAYG pensions. We find that steady-state welfare effects are 38% of consumption. We also compute the transition that follows an unexpected abolition of the PAYG pension system. This amounts to the solution of a non-linear first-order difference equation in two variables, the capital stock k_t and labor l_t. The boundary conditions are provided by the initial and final values of the capital stock, k_0 and k_T. In order to solve this numerical problem, we have to use the method of REVERSE SHOOTING described in Appendix 4.1 of my book on 'Public Economics' (2019).
We also have to use a nested version of the non-linear equations solver 'fsolve'. In particular, in the outer loop we compute the value of the capital stock in the-second-to-last period, k_T-1, that implies a value of k_0 (by backward iteration) equal to the initial steady state value of the capital stock. Within this function 'findkt()', we solve another set of non-linear equations 'kdyn()' in each period t during the backwards iteration that solves for k_t-1 and l_t-1 given k_t and l_t.
The Python command 'fsolve' is very convenient. It allows us to economize on the use of globals in the functions, which you should always avoid to comply with good programming practices. The command 'fsolve(f,x,args(y))' allows you to find the root of a function f(x,y) with respect to the variable x. In addition, you can pass on the global variables y to this function and avoid that the global variables y get overwritten.
The following Python code is stored in the file 'Ch6_social_security2.py' and can be downloaded from my web page as well:
https://www.uni-augsburg.de/de/fakultaet/wiwi/prof/vwl/heer/pubec-buch/computer-code/
# part 1: import libraries
import scipy.optimize
import numpy as np
import matplotlib.pyplot as plt
import math
plt.rcParams['figure.figsize'] = (10,6)
log = math.log
e = math.e
# part 2: parameterization
lss = 0.30 # labor supply in initial steady state (calibration of nu0)
alpha = 0.36 # production elasticity of capital
b = 0.40 # discount factor
n = 0.10 # population growth rate
nu1 = 3.33 # inverse of Frisch elasticity
tau = (0, 0.30) # parameter values for pension replacement rate (tuple)
nt=5 # number of transition periods
# part 3: defining of the functions
#
# wage function
def wage(k,l):
return (1-alpha) * k**alpha * l**(-alpha)
# interest rate function
def interest_rate(k,l):
return alpha * k**(alpha - 1) * l**(1-alpha)
# steady state computation
# input: guess for k
# output: steady-state condition (=0 in steady state)
def ksteady(x):
k = x[0]
labor = x[1]
# initialization of the return function
y = np.zeros(2)
w = wage(k,labor) # wage
r = interest_rate(k,labor) # interest rate
d = tau0 * w * labor # pension
# non-linear equations (6.28) and (6.30) in Heer (2019)
y[0] = (1+n)*k - (1-tau0)*w*labor*b/(1+b) + 1/(1+b)*(1+n)/(1+r)*d
y[1] = (1-tau0)*w - nu0*labor**nu1 /(1+b)*((1-tau0)*w*labor+(1-n)/(1+r)*d)
return y
# kdyn(x)
# procedure that computes [k_t+1,l_t+1] given [k_t,l_t]
# solution to the non-linear eqs (6.28) and (6.30) in Heer (2019)
# input: guess x0 = [k_t+1,l_t+1]
# output: the value of the equilibrium condition
# global variables: k0,l0 with are the values for k_t, l_t in the eqs
def kdyn(x0,k1,l1,tau0,tau1):
y = np.zeros(2)
k0 = x0[0]
l0 = x0[1]
r1 = interest_rate(k1,l1)
w0 = wage(k0,l0)
w1 = wage(k1,l1)
d1 = tau1*w1*l1
y[0] = k1*(1+n)-b/(1+b)*(1-tau0)*w0*l0+1/(1+b)*(1+n)/(1+r1)*d1
y[1] = (1-tau0)*w0-nu0*l0**nu1/(1+b)*( (1-tau0)*w0*l0 + (1+n)/(1+r1)*d1 )
return y
# findl(x)
# computes the optimal labor supply l_t from the foc (6.28)
def findl(l0,k0,k1,l1,tau0,tau1):
w0 = wage(k0,l0)
w1 = wage(k1,l1)
r1 = interest_rate(k1, l1)
d1 = tau1 * w1 * l1
y = (1-tau0)*w0 - nu0*l0**nu1 / (1+b)*( (1-tau0)*w0*l0 + (1+n)/(1+r1)*d1 )
return y
# findkt(x)
# computes k_0 by iterating backwards over (6.28) and (6.30)
# given a guess for k_T in the final period of the transition
# and compares the imputed value with the value in the initial
# steady state
# input: k_T
# output: k_0 - kss
# global variables that are overwritten: kt, labort, index, k0, l0
#
def findkt(x):
# value of k in the last period; afterwards k is equal to kss
kt[nt] = x
# we need to find labort[nt]
# the optimal labor supply follows from the first-order condition of labor
tau0 = tau[0]
tau1 = tau[0]
k0 = x
k1 = kss
l1 = lss
l0 = scipy.optimize.fsolve(findl, lss, args=(k0,k1,l1,tau0,tau1))
labort[nt] = l0 # labor supply in last period of transition
for i in range(nt,0,-1):
# i=0: Period 0 with steady state, tau=30% for young and old
# i=1: Period 1, young generation still has to finance old agents
# tau0=30% so that d_1=tau_1 w_1 l, but will not receive
# a pension in old age in period 2 (i=3) with tau1=0% so that d_2=0
if i==1:
tau1 = tau[1]
if i==2:
tau0 = tau[1]
k1 = kt[i]
l1 = labort[i]
x0 = [k1,l1]
y0 = scipy.optimize.fsolve(kdyn, [k1,l1], args=(k1,l1,tau0,tau1))
kt[i-1] = y0[0]
labort[i-1] = y0[1]
return kt[0]-kssp
# part 4: main program
# 4.1 Calibration of the initial steady state: parameter nu0
tau0 = tau[0] # initial steady state: no pension
kss = ( b/(1+b) * (1-alpha) / (1+n) )**(1/(1-alpha)) * lss
wss = wage(kss,lss)
rss = interest_rate(kss,lss)
c1ss = 1/(1+b)*wss*lss
nu0 = wss / ( lss**nu1 * c1ss)
c2ss = b*(1+rss)*c1ss
utilss = log(c1ss) + b*log(c2ss) - nu0*lss**(1+nu1) / (1+nu1)
print("Calibration of nu0: " + str(nu0))
# test if steady-conditions are specified correctly (=0)
x = [kss, lss]
print("steady state without pensions:")
print(x)
# 4.2 computation of the final steady state with PAYG pension
tau0 = tau[1]
# steady state with pensions
x = scipy.optimize.fsolve(ksteady, x)
kssp = x[0]
lssp = x[1]
print("solution k and l in steady state with pensions: ")
print(x)
# new steady state values of w, r, c1 and c2
wssp = wage(kssp,lssp)
rssp = interest_rate(kssp,lssp)
c1ssp = 1/(1+b) * ( (1-tau0)*wssp*lssp+(1+n)/(1+rssp)*tau0*wssp*lssp )
c2ssp = b * (1+rssp) * c1ssp
utilssp = log(c1ssp) + b * log(c2ssp) - nu0*lssp**(1+nu1) / (1+nu1)
# consumption equivalent change in steady state
cec = e**( (utilssp-utilss)/(1+b) ) - 1
print("welfare loss: " +str(100*cec) + " % of consumption")
# 4.3 computation of the dynamics: Reverse Shooting
kt = np.zeros(nt+1)
labort = np.zeros(nt+1)
periods = np.linspace(0,nt,nt+1)
# initial guess for k_nt in final period of transition
kinit = 0.999*kss
ksolution = scipy.optimize.fsolve(findkt, kinit)
print("solution found? findkt(ksolution)==0?")
print(findkt(ksolution))
fig, axes = plt.subplots(2, 1, figsize=(8, 16))
axes[0].set_xlabel('time')
axes[0].set_ylabel('capital')
axes[0].plot(kt)
axes[1].set_xlabel('time')
axes[1].set_ylabel('labor')
axes[1].plot(labort)
plt.show()
Calibration of nu0: 257.15383339746097 steady state without pensions: [0.01817586440831432, 0.3] solution k and l in steady state with pensions: [0.00674905 0.29493248] welfare loss: -38.46581175728853 % of consumption solution found? findkt(ksolution)==0? 4.8013676368086067e-14
The output of the Python code provides the steady-state values of the capital stock and labor and the welfare effects (equal to 38% of total consumption) associated with a long-run abolition of the public pension system.
The two figures illustrate the transition dynamics of labor and capital if pensions are abolished in period 1 (unexpectedly).
If you enjoyed this introductory tutorial on Python code and would like to find out more about the economics you might find it helpful to watch my youtube graduate course on 'Social Security in Dynamic General Equilibrium':
https://www.youtube.com/watch?v=0Lu2KAWLQHE&list=PL-anmMgSYtuFJu4A1biSqpKk-Oh5woKQn
Last update: Dec 11, 2020