-
Notifications
You must be signed in to change notification settings - Fork 131
/
RBC_Hansl.inp
133 lines (96 loc) · 4.45 KB
/
RBC_Hansl.inp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
## Basic RBC model with full depreciation
#
# Federico Giri and Riccardo (Jack) Lucchetti
# Ancona, 2019-01-09
#
# (adapted for hansl from JFV's matlab file)
# --- 0. Housekeeping ------------------------------------------------------
clear
set echo off
set messages off
set stopwatch
# --- 1. Calibration -------------------------------------------------------
scalar aalpha = 1/3 # Elasticity of output w.r.t. capital
scalar bbeta = 0.95 # Discount factor
# Productivity values
matrix vProductivity = {0.9792, 0.9896, 1.0000, 1.0106, 1.0212}
# Transition matrix
matrix mTransition = {0.9727, 0.0273, 0.0000, 0.0000, 0.0000; \
0.0041, 0.9806, 0.0153, 0.0000, 0.0000; \
0.0000, 0.0082, 0.9837, 0.0082, 0.0000; \
0.0000, 0.0000, 0.0153, 0.9806, 0.0041; \
0.0000, 0.0000, 0.0000, 0.0273, 0.9727}
# --- 2. Steady State ------------------------------------------------------
scalar capitalSteadyState = (aalpha*bbeta)^(1/(1-aalpha))
scalar outputSteadyState = capitalSteadyState^aalpha
scalar consumptionSteadyState = outputSteadyState-capitalSteadyState
printf " Output = %2.6f, Capital = %2.6f, Consumption = %2.6f\n", \
outputSteadyState, capitalSteadyState, consumptionSteadyState
printf "\n"
# We generate the grid of capital
scalar kmin = 0.5*capitalSteadyState
scalar kmax = 1.5*capitalSteadyState
scalar Kgrid_eps = 0.00001
scalar nGridCapital = floor((kmax - kmin)/Kgrid_eps) + 1
vGridCapital = kmin + seq(0, nGridCapital-1) * Kgrid_eps
printf "min(K) = %16.13f, max(K) = %16.13f, ngrid = %d\n", \
vGridCapital[1], vGridCapital[nGridCapital], nGridCapital
nGridProductivity = cols(vProductivity)
# --- 3. Required matrices and vectors ------------------------------------
matrix mOutput = zeros(nGridCapital,nGridProductivity)
matrix mValueFunction = zeros(nGridCapital,nGridProductivity)
matrix mValueFunctionNew = zeros(nGridCapital,nGridProductivity)
matrix mPolicyFunction = zeros(nGridCapital,nGridProductivity)
matrix expectedValueFunction = zeros(nGridCapital,nGridProductivity)
# --- 4. We pre-build output for each point in the grid ------------------
matrix mOutput = (vGridCapital'.^aalpha)*vProductivity
# --- 5. Main iteration --------------------------------------------------
scalar maxDifference = 10.0
scalar tolerance = 0.0000001
scalar iteration = 0
loop while maxDifference > tolerance --quiet
expectedValueFunction = mValueFunction*mTransition'
loop p = 1 .. nGridProductivity --quiet
# We start from previous choice (monotonicity of policy function)
scalar gridCapitalNextPeriod = 1
loop k = 1 .. nGridCapital --quiet
scalar valueHighSoFar = -1000.0
capitalChoice = vGridCapital[1]
loop n = gridCapitalNextPeriod .. nGridCapital --quiet
scalar consumption = mOutput[k,p] - vGridCapital[n]
scalar valueProvisional = (1-bbeta) * log(consumption) + bbeta * expectedValueFunction[n,p]
if valueProvisional > valueHighSoFar
valueHighSoFar = valueProvisional
capitalChoice = vGridCapital[n]
gridCapitalNextPeriod = n++
else
break # We break when we have achieved the max
endif
endloop
mValueFunctionNew[k,p] = valueHighSoFar
mPolicyFunction[k,p] = capitalChoice
endloop
endloop
maxDifference = maxr(maxc(abs(mValueFunctionNew-mValueFunction)))
mValueFunction = mValueFunctionNew
mValueFunctionNew = zeros(nGridCapital,nGridProductivity)
iteration++
if (iteration%10)==0 || iteration == 1
printf " Iteration = %3d, Sup Diff = %2.8f\n", iteration, maxDifference
flush
endif
endloop
printf " Iteration = %d, Sup Diff = %2.8f\n", iteration, maxDifference
printf "\n"
printf " My check = %2.6f\n", mPolicyFunction[1000,3]
printf "\n"
printf "Elapsed time = %g seconds\n", $stopwatch
# --- 6. Plotting results -------------------------------------------------
matrix X = mValueFunction ~ vGridCapital'
gnuplot 1 2 3 4 5 6 --matrix=X --with-lines --output=display
matrix X = mPolicyFunction ~ vGridCapital'
gnuplot 1 2 3 4 5 6 --matrix=X --with-lines --output=display
matrix vExactPolicyFunction = aalpha*bbeta.*(vGridCapital.^aalpha)
matrix dif = vExactPolicyFunction' - mPolicyFunction[,3]
matrix X = 100* dif ./ mPolicyFunction[,3]
gnuplot 1 --matrix=X --time --with-lines --output=display