Use Python to develop DSGE model

1352

收藏 2021-04-04

<div class="pinggu_markdown">
<div class="pinggu_markdown__html"><blockquote>
do it best , economy and management. 
Author : Daniel tulips liu . 
University : Renmin Universtiy of China. 
Copyright © heerburk
</blockquote>
<h2 id="a-dsge-model-use-python">A DSGE Model use Python</h2>
我了显示的格式更美观，必须要打字一部分，占用帖子的前面部分。 
很久没有能发 MARKDOWN 格式的帖子了， 
没想到今天又能发；
挺不错的。
嘿嘿。
OLG DSGE Model
<h2 id="codes">codes</h2>
<pre class=" language-c"><code class="prism language-c"># -*- coding: utf-8 -*-
"""
Created on Fri Dec 11 11:16:54 2020

@author: heerburk

OLG_simple_VI.py
simple value function iteration for the optimization paper with Alfred Maussner

you may need to adjust the load path for the excel file below (line 20)

"""

# Part 1: import libraries 
import pandas as pd
import numpy as np
import quantecon as qe
from scipy.stats import norm
import time
import math

# abbreviations 
exp = np.e
log = math.log

# Part 2: Import data 
data = pd.read_excel (r'C:\Users\Administrator\Desktop\upen_jesusfv 2021\python_chapter10_heer_maussner_dsge\survival_probs.xlsx')
df = pd.Datafr ame(data, columns= ['sp1','ef'])
#print(df) 
arr = np.array(df)
#print(arr) 
sp1 = arr[:,0]
ef = arr[0:45,1]
#print(ef) 

# Part 3: Parameterization of the model 
# demographics 
nage=70 # maximum age nw=45 # number of working years Rage=46 # first period of retirement nr=nage-Rage+1 # number of retirement years
popgrowth0 = 0.0075400000 # preferences
b=1.011 # discount factor sigma=2.0 # coefficient of relative risk aversion hbar=0.3 # fixed shift length mu0=123 # disutility of labor
mu1=0.5 # Frisch elasticity of labor supply
averagehours=0.3 # average number of working hours to calibrate mu0

# production 
r=0.04 # initial value of the interest rate (in order to compute kbar)
alp=0.36 # production elasticity of capital delta=0.06 # rate of depreciation
# productivity of workers 
lamb0=.96 # autoregressive parameter sigmay1=0.38 # variance for 20-year old, log earnings */
sigmae=0.045 # earnings disturbance term variance */
ny=5 # number of productivity types
m=2 # width of the productivity grid

# parameters of the algorighm 
psi1=0.0001 # parameter of utility function for zero or negative consumption
update_agg=0.8 # parameter for the update of initial values
nq=50 # maximum number of iterations to find a solution

tol=0.001 # percentage deviation of final solution
# grids 
kmin=0 # inidividual wealth
kmax=20 # upper limit of capital grid na=100 # number of grid points over assets a in [kmin,kmax]
a = np.linspace(kmin, kmax, na) # asset grid emin=0 # individual cumulated earnings
emax=3
nce=20 # number of grid points on cumulated earnings
e = np.linspace(emin, emax, nce) # grid over cumulated earnings
labormin=0
labormax=0.7 # maximum labor supply
nlabor=100
lgrid= np.linspace(labormin, labormax, nlabor) # grid over labor supply
ymin=0
ymax=4
nin=100
incomegrid = np.linspace(ymin, ymax, nin) # grid for gross income

# measure of living persons 
mass = np.ones(nage)

for i in range(nage-1):
print(i)
mass[i+1]=mass[i]*sp1[i]/(1+popgrowth0)

mass = mass / mass.sum()

# compute productivity grid 
sy = np.sqrt(sigmae/(1-lamb0**2)) # earnings variance x0=1 nsim=10
ye = np.linspace(-m*sy, m*sy, ny)
ye1 = exp**ye

# transition matrix using Tauchen's approximation 
# return is a class 'Markov Chain' 
mc = qe.markov.approximation.tauchen(lamb0, sigmae, 0.0, m, ny)
# transition matrix is stored in ob ject 'PÄ 
py = mc.P

# mass of the workers 
muy = np.zeros((nw,ny))
w = ye[1]-ye[0]

# first year mass distribution 
muy[0,0] = norm.cdf( (ye[0]+w/2)/np.sqrt(sigmay1) ) * mass[0]
muy[0,ny-1] = 1-norm.cdf( (ye[ny-1]-w/2)/np.sqrt(sigmay1))
muy[0,ny-1] = muy[0,ny-1]*mass[0]

for i in range(1,ny-1):
print(i)
muy[0,i] = ( norm.cdf( (ye[i]+w/2)/np.sqrt(sigmay1) ) -
 norm.cdf( (ye[i]-w/2)/np.sqrt(sigmay1) ) )
muy[0,i] = muy[0,i]*mass[0]


for i in range (1,nw):
print(i)
muy[i,:] = muy[i-1,:] @ py * sp1[i] / (1+popgrowth0)

# initial guess for aggregate effective labor supply 
temp = np.transpose(muy*ye1) * ef
nbar = temp.sum() * hbar / sum(mass[0:nw])

xmat = np.zeros((4,nq)) # stores initial values at each iteration q

# initial guesses 
kbar = (alp/(r+delta))**(1/(1-alp))*nbar
wage = (1-alp)*kbar**alp * nbar**(-alp)
tau = 0.068/(1-alp) # pensions in the US are 6.8 of GDP in 2012, since wage income is (1-alpha)*GDP, 6.8%/(1-alp) provides an initial guess
tau = 0.124 # social security tax in the US 2020 trbar = np.mean(1-sp1[0:nage-1]) * kbar # proxy of accidental bequests

 # part 4: definition of the functions 
#
# wage function 
def wage_rate(k,l):
return (1-alp) * k**alp * l**(-alp)

# interest rate function 
def interest_rate(k,l):
return alp * k**(alp - 1) * l**(1-alp) - delta


# utility function # I changed this: only if c.le.0, a small amout is added @ 
def u(x,y):
x1 = np.where(x<=0, psi1, x) # replace negative values by psi1
if sigma==1:
 return log(x1) - mu0*y**(1+1/mu1) / (1+1/mu1)
else:
 return x1**(1-sigma) / (1-sigma) -mu0*y**(1+1/mu1) / (1+1/mu1)

# pension schedule 
# see Huggett and Parra (2010, JPE)def pension(x): 
bendpoint1 = 0.20*meanincome
bendpoint2 = 1.24*meanincome
bmin = penmin


if x<bendpoint1:
 y = bmin+0.9*x
elif x>=bendpoint1 and x<=bendpoint2:
 y = bmin + 0.9*bendpoint1 + 0.32*(x-bendpoint1)
else:
 y = bmin+ 0.9*bendpoint1 + 0.32*(bendpoint2-bendpoint1) + 0.15*(x-bendpoint2)

return y

# computes the gini for distribution where x has measure g(x) def ginid(x,gy,ng): 
x = np.where(x<=0, 0, x)
xmean = np.dot(x,gy)
# y = np.c_[x,gy] 
# sorting according to first column still has to be implemented# y0 = np.sort(y,axis=0) 
# x = y0[:,0] 
# g = y0[:,1] 
 g = gy
f = np.zeros(ng) # accumulated frequency */
 f[0] = g[0]*x[0]/xmean
gini = 1- f[0] * g[0]
for i in range(1,ng):
 f[i] = f[i-1] + g[i]*x[i]/xmean
 gini = gini - (f[i]+f[i-1])*g[i]

return gini

start_time = time.time()

# part 5: main program 
# outer loop 
q=0
crit = 1+tol
while q<nq and crit>tol:
print("q: " + str(q))
xmat[0,q] = kbar
xmat[1,q] = nbar
xmat[2,q] = tau
xmat[3,q] = trbar


# compute variables derived from (kbar,nbar,tau,trbar) 
# required for the computation of decision rules 
wage = wage_rate(kbar, nbar)
rate = interest_rate(kbar, nbar)
ybar=kbar**(alp) * nbar**(1-alp) # produciton
 meanincome=wage*nbar / np.sum(mass[0:nw]) # average earnings of the workers
 penmin=0.1242*ybar # minimum pension
 start_time1 = time.time() # time to compute decision rules
# initialization of agents' value function: retired households 
 vr = np.zeros((na*nce,nr)) # value function
aropt = np.zeros((na*nce,nr)) # optimal asset
cropt = np.zeros((na*nce,nr)) # optimal consumption
 k0 = 0 # initialization of workers
 k1 = 0 # current assets and cumulated earnings,
 e0 = 0 # next-period assets and
e1 = 0 # cumulated earnings
</code></pre>
<h1 id="enjoy">Enjoy</h1>
Yours Daniel Tulips Liu
Renmin University of China , Beijing.
</div>
</div>

扫码加我拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

扫码加我 拉你入群

分享

扫码加好友，拉您进群

扫码加我拉你入群