-
Notifications
You must be signed in to change notification settings - Fork 2
/
optimal_path.m
150 lines (132 loc) · 4.36 KB
/
optimal_path.m
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
function results = optimal_path(parameters,f0)
% This version codes Galo's solution using standard finite difference methods.
% Transition path, no shocks
%% Model parameters
% tolerance
tol = parameters.tol;
tol_path = parameters.tol_path;
max_iter = parameters.max_iter; % Maximum number of iterations
relax = parameters.relax; % Relaxation coefficient
relax_path = parameters.relax_path;
% numerical parameters
t_max = parameters.t_max; % Maximum number of years (time)
tau_max = parameters.tau_max; % Maximum number of years (maturity)
dt = parameters.dt; % Monthly steps
% preference Parameters
rho = parameters.rho; % Discount factor - Model is Quarterly
delta = parameters.delta; % Coupon portion
gamma = parameters.gamma; % risk aversion
% frictions parameters
lambda_bar = parameters.lambda_bar;
% psi
psi_n=parameters.psi;
psi = psi_n(:,end);
%% Variables (preallocation)
t = 0:dt:t_max; % calendar time
N = length(t); % Number of time nodes
tau = dt:dt:tau_max; % maturities
I = length(tau); % Number of tau nodes
u = zeros (I,1);
u(:) = -delta;
u(1) = u(1) -1/dt;
A0 = -speye(I);
aa = ones(I-1,1);
A1 = spdiags(aa,-1,I,I);
A = 1/dt*(A0+A1);
V_n = zeros(I,N); % Value function
iota_n = zeros(I,N); % Policy function
f_n = zeros(I,N); % Density
c_n = zeros(1,N); % Consumption
r_n_u = zeros(1,N); % Interest rate (update)
%Path for Income
y_bar = parameters.y_bar ; %Before the Shock
y_bar_0 = parameters.y_bar_0;
rho_y = parameters.rho_y;
y_n = zeros(1,N);
for i=2:N
y_n(1) = y_bar_0;
y_n(i) = y_bar *(1-rho_y)+rho_y*y_n(i-1);
end
%% Initial guess of interest rates
r_n = rho * ones(1,N);
%% STEADY STATE
V = -delta/rho *(1-exp(-rho*tau'))-exp(-rho*tau'); % HJB: Solution at Steady State
iota = 1/lambda_bar * (V+ psi );% Optimal policies
f = flip(cumsum(flip(iota)))*dt;% KFE
c = y_bar - f(1) + sum((psi - 1/2 *lambda_bar * iota).*iota-delta*f) * dt; % Consumption
r = rho;
if c<=0 % The case lambda <lambda_min
r = (r_bar + rho) /2; % First guess
for iter = 1:max_iter
V = -delta/r *(1-exp(-r*tau'))-exp(-r*tau');
iota = 1./lambda_bar .* (V+ psi );% Optimal policies
f = flip(cumsum(flip(iota)))*dt;% KFE
c = y_bar - f(1) + sum((psi - 1/2 *lambda_bar .* iota).*iota-delta*f) * dt;% Consumption
if abs(c) < tol_path
break % Consumption should be zero
else
r = r + relax_path*c;
end
end
if iter == max_iter
disp('Error: maximum number of iterations reached at steady state loop')
end
end
%% DYNAMICS
%%----------------------------
% Initial guess of interest rates
r_n = r * ones(1,N);
tic
for iter = 1:max_iter
% HJB: Dynamic solution
V_n(:,N) = V; % Terminal value
for n = N-1:-1:1
B = (1/dt + r_n(n))*speye(I) -A;
d = u + V_n(:,n+1) /dt;
V_aux = B\d;
V_n(:,n) = V_aux;
end
% Optimal policies
iota_n = 1/lambda_bar * (V_n + psi_n );
% KFE
f_n(:,1) = f0;
D = speye(I) - dt *A';
for n=2:N
h = iota_n(:,n) *dt + f_n(:,n-1);
f_aux = D\h;
f_n(:,n) = f_aux;
end
% Consumption and interest rates
c_n = y_n - f_n(1,:) + sum((psi_n - 1/2 *lambda_bar * iota_n).*iota_n-delta*f_n) * dt;
if c>tol_path
r_n_u(1:N-1) = rho * ones(1,N-1) + gamma/dt * (c_n(2:N) - c_n(1:N-1)) ./ c_n(1:N-1);
r_n_u(N) = rho;
else
r_n_u(1:N-2) = rho * ones(1,N-2) + gamma/dt * (c_n(2:N-1) - c_n(1:N-2)) ./ c_n(1:N-2);
r_n_u(N-1) = r;
r_n_u(N) = r;
end
% Interest rate update
error = max(abs (r_n - r_n_u));
if error < tol_path
break
else
r_n = relax_path * r_n_u + (1-relax_path) * r_n;
if mod(iter,500) == 0
disp(iter)
disp(error)
end
end
end
if iter == max_iter
disp('Error: maximum number of iterations reached')
end
toc
results.V_n = V_n; % Value function
results.iota_n = iota_n;% Policy function
results.f_n = f_n; % Density
results.c_n = c_n; % Consumption
results.r_n = r_n; % Interest rate (update)
results.y_n = y_n;
results.psi = psi_n;
results.dt = dt;