-
Notifications
You must be signed in to change notification settings - Fork 2
/
OBC_App_2_sims_PF.m
105 lines (78 loc) · 3.13 KB
/
OBC_App_2_sims_PF.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
%Policy Application from the OBC paper (NK model), perfect forsight paths.
%To study a different example, simply change 'Insert' files and nvar,nx
%Model structures are defined in the 'Insert' files
%Written by Michael Hatcher ([email protected]). Any errors are my own.
clc; clear; %close all;
T_guess = 50; % Final date before terminal solution (guess)
T_sim = T_guess + 1; % Simulation length
T_news = 1;
N_guess = 100; %No. of guesses
%Shocks
zero = zeros(1,T_sim+1);
% Model and calibration
run Insert_App_2
%run Insert_App_2_PLT
%No. of variables
nvar = length(B1(:,1));
nx = 0; % %No. exogenous variables in x
% Find terminal solution (Cho and Moreno 2011, JEDC)
run Cho_and_Moreno
Omeg_t = NaN(size(Omega_bar,1), size(Omega_bar,2), T_sim);
Gama_t = NaN(size(Gama_bar,1), size(Gama_bar,2), T_sim);
Psi_t = NaN(size(Psi_bar,1), size(Psi_bar,2), T_sim);
%Guessed structure
rng(1)
ind_stack = randi([0 1],T_guess,N_guess);
ind_stack(:,1) = ones(T_guess,1);
for i=2:T_guess
ind_stack(:,i) = [zeros(i-1,1); ones(T_guess-(i-1),1)];
end
vec_1 = ones(T_sim-T_guess,1);
%Check if M is a P matrix
e(1) = 0.01; e(2:T_news) = 0; %Specified news shocks
e_vec = [e(1) e(2:T_news) zeros(1,T_sim+1-T_news)];
X_init = zeros(nvar,1); %Initial values
%run M_matrix
%run P_matrix
mstar = zeros(N_guess,1);
run PF_insert.m
mstar(mstar==0) = [];
if isempty(mstar)
disp('No solution found. Check T_guess and N_guess are not too small.')
elseif ~isempty(mstar)
X_sol = X_sol(:,:,mstar);
X_exog = X_sol(nvar-nx+1:nvar,:,1);
X_sol_exc = X_sol_exc(:,:,mstar);
%Note: The below is an ad-hoc workaround. Because in the above model i(t) = i*(t) in the reference regime,
%ones of these paths will removed by the 'unique' function. To prevent this we add small positive number to
%final simulated point of variable 1 (i.e. i(t)).
for ii=1:length(X_sol(1,1,:))
orig = X_sol(1,T_sim-1,ii); orig_exc = X_sol_exc(1,T_sim-1,ii);
X_sol(1,T_sim-1,ii) = orig + (1+abs(rand)/100)*1e-15;
X_sol_exc(1,T_sim-1,ii) = orig_exc + (1+abs(rand)/100)*1e-15;
end
solutions = reshape(permute(X_sol_exc,[1,3,2]),[],size(X_sol_exc,2));
ind_solutions = ind_sol(mstar,:);
%Solution paths
x_fin = unique(solutions,'rows','stable');
x_fin = [x_fin; X_exog];
sol_ind = unique(ind_solutions,'rows','stable');
ind_fin = sol_ind(1,:);
%No. of solutions
k = length(x_fin(:,1))/nvar;
if k==1
disp('Unique solution found')
else
disp('Multiple solutions found')
end
end
figure(1)
T_plot = 31; %Time horizon in plot
T_sim2 = (T_sim<=T_plot)*T_sim + (T_sim>T_plot)*T_plot;
Time = 1:T_sim2-1; %Counter for time
for j=1:k
subplot(1,3,1), plot(Time,x_fin(4+nvar*(j-1),1:T_sim2-1),'--g'), xlabel('Time'), title('Inflation'), xlim([1 inf]), hold on
subplot(1,3,2), plot(Time,x_fin(3+nvar*(j-1),1:T_sim2-1),'--g'), xlabel('Time'), title('Output gap'), xlim([1 inf]), hold on
subplot(1,3,3), plot(Time,100*(x_fin(1+nvar*(j-1),1:T_sim2-1)-X1_min),'--g'), xlabel('Time'), title('Interest rate (%)'), xlim([1 inf]), hold on,
end
ylim([-inf inf])