-
Notifications
You must be signed in to change notification settings - Fork 0
/
Mathematica_Plain_Text.txt
127 lines (100 loc) · 3.71 KB
/
Mathematica_Plain_Text.txt
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
Clear["Global`*"];
AbsoluteTiming[
(* Compiling the inner Loop *)
(* NOTE: Using the un-documented Compile`GetElement function instead \
of Part *)
innerLoop =
Compile[{{mOutput, _Real, 2}, {vGridCapital, _Real,
1}, {nGridCapital, _Integer},
{nGridProductivity, _Integer}, {expectedValueFunction, _Real,
2}},
With[{initialCapital = First[vGridCapital]},
Table[
Module[{gridCapitalNextPeriod = 1},
Table[
With[{output =
Compile`GetElement[mOutput, nCapital, nProductivity]},
Module[{
valueProvisional = 0.,
valueHighSoFar = -1000.0,
capitalChoice = initialCapital},
Catch@Do[
valueProvisional =
0.05 Log[
output -
Compile`GetElement[vGridCapital, nCapitalNextPeriod]] +
0.95 Compile`GetElement[expectedValueFunction,
nCapitalNextPeriod, nProductivity];
If[valueProvisional > valueHighSoFar,
valueHighSoFar = valueProvisional;
capitalChoice =
Compile`GetElement[vGridCapital, nCapitalNextPeriod];
gridCapitalNextPeriod = nCapitalNextPeriod;
,
Throw[{valueHighSoFar, capitalChoice}]
],
{nCapitalNextPeriod, gridCapitalNextPeriod,
nGridCapital}]]],
{nCapital, nGridCapital}]],
{nProductivity, nGridProductivity}]
],
CompilationTarget -> "C", "RuntimeOptions" -> "Speed" ];
(* 1. Calibration*)
\[Alpha] = 0.333333333333;
\[Beta] = 0.95;
(* Productivity values*)
vProductivity = {0.9792, 0.9896, 1.0000, 1.0106, 1.0212};
(* Transition 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}};
mTransitionTransposed = Transpose[mTransition];
(* 2. Steady State*)
Subscript[k, ss] = (\[Alpha] \[Beta])^(1/(1 - \[Alpha]));
Subscript[y, ss] = Subscript[k, ss]^\[Alpha];
Subscript[c, ss] = Subscript[y, ss] - Subscript[k, ss];
(* We generate the grid of capital*)
vGridCapital =
Range[0.5 Subscript[k, ss], 1.5 Subscript[k, ss], 0.00001];
nGridCapital = Length[vGridCapital];
nGridProductivity = Length[vProductivity];
(*3. Required matrices and vectors*)
mOutput = ConstantArray[0, {nGridCapital, nGridProductivity}];
mValueFunction =
ConstantArray[0, {nGridCapital, nGridProductivity}];
mValueFunctionNew =
ConstantArray[0, {nGridCapital, nGridProductivity}];
mPolicyFunction =
ConstantArray[0, {nGridCapital, nGridProductivity}];
expectedValueFunction =
ConstantArray[0, {nGridCapital, nGridProductivity}];
(*4. We pre-build output for each point in the grid*)
mOutput = Transpose[{vGridCapital^\[Alpha]}].{vProductivity};
(* FixedPoint *)
tolerance = 0.0000001;
iteration = 0;
dis = 0;
(* outer Loop function *)
outerLoop[{mValueFunction_, mPolicyFunction_}] := Transpose[
innerLoop[mOutput, vGridCapital, nGridCapital, nGridProductivity,
Dot[mValueFunction, mTransitionTransposed]],
{3, 2, 1}];
(* Iteration *)
{mValueFunction, mPolicyFunction} =
FixedPoint[outerLoop, {mValueFunction, mPolicyFunction},
SameTest ->
(
(dis = Max[Abs[#1[[1]] - #2[[1]]]];
iteration++;
If[Mod[iteration, 10] == 0 || iteration == 1,
Print[
StringForm["Iteration = ``, Sup Diff = ``", iteration,
dis]]];
dis < tolerance ) &
)
];
Print[StringForm["Iteration = ``, Sup Diff = ``", iteration, dis]];
Print[StringForm["My check = ``", mPolicyFunction[[1000, 3]]]];
]