diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 5af50faa2733af6b162b2bdf53e44d596f64a64a..66871d8dc66e20fe594bb2a0cbb9ec1895c3b3f8 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -357,7 +357,6 @@ options_.drop = 100;
options_.aim_solver = false; % i.e. by default do not use G.Anderson's AIM solver, use mjdgges instead
options_.k_order_solver = false; % by default do not use k_order_perturbation but mjdgges
options_.partial_information = false;
-options_.ACES_solver = false;
options_.conditional_variance_decomposition = [];
% Ramsey policy
diff --git a/matlab/partial_information/PCL_Part_info_moments.m b/matlab/partial_information/PCL_Part_info_moments.m
index 60f4a9bc9a50f68adb07ceb3f8e1b01f02ea3384..f4932d3f5efd01bd7098d0949b30de59e78640c8 100644
--- a/matlab/partial_information/PCL_Part_info_moments.m
+++ b/matlab/partial_information/PCL_Part_info_moments.m
@@ -1,11 +1,12 @@
-function AutoCOR_YRk=PCL_Part_info_moments( H, varobs, dr,ivar)
+function oo_=PCL_Part_info_moments(M_, oo_, options_, varobs, dr, ivar)
+% function oo_=PCL_Part_info_moments(M_, oo_, options_, varobs, dr, ivar)
% sets up parameters and calls part-info kalman filter
% developed by G Perendia, July 2006 for implementation from notes by Prof. Joe Pearlman to
% suit partial information RE solution in accordance with, and based on, the
% Pearlman, Currie and Levine 1986 solution.
% 22/10/06 - Version 2 for new Riccati with 4 params instead 5
-% Copyright © 2006-2018 Dynare Team
+% Copyright © 2006-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -27,7 +28,10 @@ function AutoCOR_YRk=PCL_Part_info_moments( H, varobs, dr,ivar)
% and the jump variables x(t).
% The jump variables have dimension NETA
-global M_ options_ oo_
+if ~exist([M_.dname '/Output'],'dir')
+ mkdir(M_.dname,'Output');
+end
+
warning_old_state = warning;
warning off
@@ -40,29 +44,16 @@ NOBS = length(OBS);
G1=dr.PI_ghx;
impact=dr.PI_ghu;
nmat=dr.PI_nmat;
-CC=dr.PI_CC;
NX=M_.exo_nbr; % no of exogenous varexo shock variables.
FL_RANK=dr.PI_FL_RANK;
NY=M_.endo_nbr;
LL = sparse(1:NOBS,OBS,ones(NOBS,1),NY,NY);
-if exist( 'irfpers')==1
- if ~isempty(irfpers)
- if irfpers<=0, irfpers=20, end
- else
- irfpers=20;
- end
-else
- irfpers=20;
-end
-
ss=size(G1,1);
pd=ss-size(nmat,1);
SDX=M_.Sigma_e^0.5; % =SD,not V-COV, of Exog shocks or M_.Sigma_e^0.5 num_exog x num_exog matrix
-if isempty(H)
- H=M_.H;
-end
+H=M_.H;
VV=H; % V-COV of observation errors.
MM=impact*SDX; % R*(Q^0.5) in standard KF notation
% observation vector indices
@@ -93,17 +84,14 @@ A11=G1(1:pd,1:pd);
A22=G1(pd+1:end, pd+1:end);
A12=G1(1:pd, pd+1:end);
A21=G1(pd+1:end,1:pd);
-Lambda= nmat*A12+A22;
-I_L=inv(Lambda);
-BB=A12*inv(A22);
-FF=K2*inv(A22);
+BB=A12/A22;
+FF=K2/A22;
QQ=BB*U22*BB' + U11;
UFT=U22*FF';
% kf_param structure:
AA=A11-BB*A21;
CCCC=A11-A12*nmat; % F in new notation
DD=K1-FF*A21; % H in new notation
-EE=K1-K2*nmat;
RR=FF*UFT+VV;
if ~any(RR)
% if zero add some dummy measurement err. variance-covariances
@@ -113,30 +101,23 @@ if ~any(RR)
RR=eye(size(RR,1))*1.0e-6;
end
SS=BB*UFT;
-VKLUFT=VV+K2*I_L*UFT;
-ALUFT=A12*I_L*UFT;
-FULKV=FF*U22*I_L'*K2'+VV;
-FUBT=FF*U22*BB';
-nmat=nmat;
% initialise pshat
AQDS=AA*QQ*DD'+SS;
DQDR=DD*QQ*DD'+RR;
-I_DQDR=inv(DQDR);
-AQDQ=AQDS*I_DQDR;
+AQDQ=AQDS/DQDR;
ff=AA-AQDQ*DD;
hh=AA*QQ*AA'-AQDQ*AQDS';%*(DD*QQ*AA'+SS');
rr=DD*QQ*DD'+RR;
ZSIG0=disc_riccati_fast(ff,DD,rr,hh);
PP=ZSIG0 +QQ;
-exo_names = M_.exo_names;
-
DPDR=DD*PP*DD'+RR;
-I_DPDR=inv(DPDR);
-PDIDPDRD=PP*DD'*I_DPDR*DD;
-MSIG=disclyap_fast(CCCC, CCCC*PDIDPDRD*PP*CCCC', options_.lyapunov_doubling_tol);
+PDIDPDRD=PP*DD'/DPDR*DD;
+
+MSIG=lyapunov_solver(CCCC,CCCC,PDIDPDRD*PP,options_);
COV_P=[ PP, PP; PP, PP+MSIG]; % P0
+COV_P=(COV_P+COV_P')/2; %assure symmetry
dr.PI_GG=[CCCC (AA-CCCC)*(eye(ss-FL_RANK)-PDIDPDRD); zeros(ss-FL_RANK) AA*(eye(ss-FL_RANK)-PDIDPDRD)];
@@ -144,10 +125,11 @@ GAM= [ AA*(eye(ss-FL_RANK)-PDIDPDRD) zeros(ss-FL_RANK); (AA-CCCC)*(eye(ss-FL_RAN
VV = [ dr.PI_TT1 dr.PI_TT2];
nn=size(VV,1);
-COV_OMEGA= COV_P( end-nn+1:end, end-nn+1:end);
+COV_OMEGA= COV_P(end-nn+1:end,end-nn+1:end);
COV_YR0= VV*COV_OMEGA*VV';
diagCovYR0=diag(COV_YR0);
-labels = M_.endo_names(ivar);
+diagCovYR0(abs(diagCovYR0)<1e-10)=0;
+labels=get_labels_transformed_vars(M_.endo_names,ivar,options_,false);
if ~options_.nomoments
z = [ sqrt(diagCovYR0(ivar)) diagCovYR0(ivar) ];
@@ -157,18 +139,36 @@ if ~options_.nomoments
end
if ~options_.nocorr
diagSqrtCovYR0 = sqrt(diagCovYR0);
- DELTA = inv(diag(diagSqrtCovYR0));
+ non_zero_indices=find(diagSqrtCovYR0>1e-10);
+ temp=1./diagSqrtCovYR0(non_zero_indices);
+ DELTA = zeros(length(diagCovYR0),1);
+ DELTA(non_zero_indices)=temp;
+ DELTA=diag(DELTA);
COR_Y = DELTA*COV_YR0*DELTA;
- title = 'MATRIX OF CORRELATION';
- headers = vertcat('VARIABLE', M_.endo_names(ivar));
- dyntable(options_, title, headers, labels, COR_Y(ivar,ivar), size(labels,2)+2, 8, 4);
+ i1=intersect(non_zero_indices,ivar);
+ if options_.contemporaneous_correlation
+ oo_.contemporaneous_correlation = COR_Y(ivar,ivar);
+ end
+ if ~options_.noprint
+ skipline()
+ title = 'MATRIX OF CORRELATIONS';
+ labels=get_labels_transformed_vars(M_.endo_names,i1,options_,false);
+ headers = vertcat('Variables', labels);
+ lh = cellofchararraymaxlength(labels)+2;
+ dyntable(options_, title, headers, labels, COR_Y(i1,i1), lh, 8, 4);
+ if options_.TeX
+ labels=get_labels_transformed_vars(M_.endo_names_tex,i1,options_,true);
+ headers = vertcat('Variables', labels);
+ lh = cellofchararraymaxlength(labels)+2;
+ dyn_latex_table(M_, options_, title, 'th_corr_matrix', headers, labels, COR_Y(i1,i1), lh, 8, 4);
+ end
+ end
else
COR_Y=[];
end
ar = options_.ar;
if ar > 0
- COV_YRk = zeros(nn, ar);
AutoCOR_YRk= zeros(nn, ar);
for k = 1:ar
COV_P = GAM*COV_P;
@@ -179,10 +179,18 @@ if ar > 0
AutoCOR_YRk(:,k) = diag(COV_YRk)./diagCovYR0;
end
title = 'COEFFICIENTS OF AUTOCORRELATION';
- headers = vertcat('VARIABLE', cellstr(int2str([1:ar]')));
- dyntable(options_, title, headers, labels, AutoCOR_YRk(ivar,:), size(labels,2)+2, 8, 4);
+ labels=get_labels_transformed_vars(M_.endo_names,i1,options_,false);
+ headers = vertcat('Order ', cellstr(int2str((1:options_.ar)')));
+ lh = cellofchararraymaxlength(labels)+2;
+ dyntable(options_, title, headers, labels, AutoCOR_YRk(i1,:), lh, 8, 4);
+ if options_.TeX
+ labels=get_labels_transformed_vars(M_.endo_names_tex,i1,options_,true);
+ headers = vertcat('Order ', cellstr(int2str((1:options_.ar)')));
+ lh = cellofchararraymaxlength(labels)+2;
+ dyn_latex_table(M_, options_, title, 'th_autocorr_matrix', headers, labels, AutoCOR_YRk(i1,:), lh, 8, 4);
+ end
else
AutoCOR_YRk = [];
end
-save ([M_.fname '_PCL_moments'], 'COV_YR0','AutoCOR_YRk', 'COR_Y');
+save ([M_.dname filesep 'Output' filesep M_.fname '_PCL_moments'], 'COV_YR0','AutoCOR_YRk', 'COR_Y');
warning(warning_old_state);
diff --git a/matlab/partial_information/PCL_resol.m b/matlab/partial_information/PCL_resol.m
index f01d30ed5eda94803b116a163c808137bd1e66aa..3ce946e47a7830503c2917ba5763b038fdde0e80 100644
--- a/matlab/partial_information/PCL_resol.m
+++ b/matlab/partial_information/PCL_resol.m
@@ -1,5 +1,5 @@
-function [dr,info]=PCL_resol(ys,check_flag)
-% function [dr,info]=PCL_resol(ys,check_flag)
+function [dr,info]=PCL_resol(M_, options_, oo_)
+% function [dr,info]=PCL_resol(M_, options_, oo_)
% Computes first and second order approximations
%
% INPUTS
@@ -10,17 +10,6 @@ function [dr,info]=PCL_resol(ys,check_flag)
% OUTPUTS
% dr: structure of decision rules for stochastic simulations
% info=1: the model doesn't determine the current variables '...' uniquely
-% info=2: MJDGGES returns the following error code'
-% info=3: Blanchard Kahn conditions are not satisfied: no stable '...' equilibrium
-% info=4: Blanchard Kahn conditions are not satisfied:'...' indeterminacy
-% info=5: Blanchard Kahn conditions are not satisfied:'...' indeterminacy due to rank failure
-% info=6: The jacobian evaluated at the steady state is complex.
-% info=19: The steadystate file did not compute the steady state (inconsistent deep parameters).
-% info=20: can't find steady state info(2) contains sum of sqare residuals
-% info=21: steady state is complex
-% info(2) contains sum of sqare of
-% imaginary part of steady state
-% info=30: Variance can't be computed
%
% SPECIAL REQUIREMENTS
% none
@@ -42,15 +31,6 @@ function [dr,info]=PCL_resol(ys,check_flag)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global M_ options_ oo_
-global it_
-
-jacobian_flag = false;
-
-info = 0;
-
-it_ = M_.maximum_lag + 1 ;
-
if M_.exo_nbr == 0
oo_.exo_steady_state = [] ;
end
@@ -62,80 +42,13 @@ if M_.exo_det_nbr > 0
tempexdet = oo_.exo_det_simul;
oo_.exo_det_simul = repmat(oo_.exo_det_steady_state',M_.maximum_lag+M_.maximum_lead+1,1);
end
-dr.ys = ys;
-check1 = 0;
-% testing for steadystate file
-static_resid = str2func(sprintf('%s.sparse.static_resid', M_.fname));
-static_g1 = str2func(sprintf('%s.sparse.static_g1', M_.fname));
-function [resid, g1] = static_resid_g1(y, x, params)
- [resid, T_order, T] = static_resid(y, x, params);
- g1 = static_g1(y, x, params, M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr, T_order, T);
-end
-
-if options_.steadystate_flag
- [dr.ys,check1] = feval([M_.fname '_steadystate'],dr.ys,...
- [oo_.exo_steady_state; ...
- oo_.exo_det_steady_state]);
- if size(dr.ys,1) < M_.endo_nbr
- if length(M_.aux_vars) > 0
- dr.ys = add_auxiliary_variables_to_steadystate(dr.ys,M_.aux_vars,...
- M_.fname,...
- oo_.exo_steady_state,...
- oo_.exo_det_steady_state,...
- M_.params);
- else
- error([M_.fname '_steadystate.m doesn''t match the model']);
- end
- end
-else
- % testing if ys isn't a steady state or if we aren't computing Ramsey policy
- if ~options_.ramsey_policy
- if options_.linear == 0
- % nonlinear models
- if max(abs(static_resid_g1(dr.ys, [oo_.exo_steady_state; ...
- oo_.exo_det_steady_state], M_.params))) > options_.solve_tolf
- opt = options_;
- opt.jacobian_flag = false;
- [dr.ys, check1] = dynare_solve(static_resid_g1, dr.ys, options.steady_.maxit, options_.solve_tolf, options_.solve_tolx, ...
- opt, [oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params);
- end
- else
- % linear models
- [fvec,jacob] = static_resid_g1(dr.ys, [oo_.exo_steady_state;...
- oo_.exo_det_steady_state], M_.params);
- if max(abs(fvec)) > 1e-12
- dr.ys = dr.ys-jacob\fvec;
- end
- end
- end
-end
-% testing for problem
-if check1
- if options_.steadystate_flag
- info(1)= 19;
- resid = check1 ;
- else
- info(1)= 20;
- resid = static_resid(ys, oo_.exo_steady_state, M_.params);
- end
- info(2) = resid'*resid ;
- return
-end
-
-if ~isreal(dr.ys)
- info(1) = 21;
- info(2) = sum(imag(ys).^2);
- dr.ys = real(dr.ys);
+[dr.ys,M_.params,info] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,~options_.steadystate.nocheck);
+if info(1)
return
end
-dr.fbias = zeros(M_.endo_nbr,1);
-if( options_.partial_information || options_.ACES_solver)%&& (check_flag == 0)
- [dr,info,M_,options_,oo_] = dr1_PI(dr,check_flag,M_,options_,oo_);
-else
- [dr,info,M_,options_,oo_] = dr1(dr,check_flag,M_,options_,oo_);
-end
+[dr,info,oo_] = dr1_PI(dr,M_,options_,oo_);
if info(1)
return
end
@@ -143,12 +56,4 @@ end
if M_.exo_det_nbr > 0
oo_.exo_det_simul = tempexdet;
end
-oo_.exo_simul = tempex;
-tempex = [];
-
-end
-
-% 01/01/2003 MJ added dr_algo == 1
-% 08/24/2001 MJ uses Schmitt-Grohe and Uribe (2001) constant correction
-% in dr.ghs2
-% 05/26/2003 MJ added temporary values for oo_.exo_simul
+oo_.exo_simul = tempex;
\ No newline at end of file
diff --git a/matlab/partial_information/PI_gensys.m b/matlab/partial_information/PI_gensys.m
index 20a860ff45dd4767a626ba7eb180888911bd4ae4..01639bf7d7076d6eeaa1da1e986b335d587802ce 100644
--- a/matlab/partial_information/PI_gensys.m
+++ b/matlab/partial_information/PI_gensys.m
@@ -1,4 +1,4 @@
-function [G1pi,C,impact,nmat,TT1,TT2,gev,eu, DD, E2, E5, GAMMA, FL_RANK ]=PI_gensys(a0,a1,a2,a3,c,PSI,NX,NETA,FL_RANK,M_,options_)
+function [G1pi,C,impact,nmat,TT1,TT2,gev,eu, DD, E2, E5, GAMMA, FL_RANK ]=PI_gensys(a0,a1,a2,c,PSI,NX)
% System given as
% a0*E_t[y(t+1])+a1*y(t)=a2*y(t-1)+c+psi*eps(t)
% with z an exogenous variable process.
@@ -35,10 +35,10 @@ function [G1pi,C,impact,nmat,TT1,TT2,gev,eu, DD, E2, E5, GAMMA, FL_RANK ]=PI_gen
lastwarn('','');
-global lq_instruments;
-eu=[0;0];C=c;
+
+C=c;
realsmall=1e-6;
-fixdiv=(nargin==6);
+fixdiv=(nargin==5);
n=size(a0,1);
DD=[];E2=[]; E5=0; GAMMA=[];
%
@@ -68,7 +68,7 @@ F3=Sinv*U01'*a2*V01;
F4=Sinv*U01'*a2*V02;
F5=Sinv*U01'*PSI;
singular=0;
-warning('', '');
+
try
if rcond(C2) < 1e-8
singular=1;
@@ -89,20 +89,12 @@ try
warning('on','MATLAB:singularMatrix');
warning('on','MATLAB:nearlySingularMatrix');
if (any(any(isinf(UAVinv))) || any(any(isnan(UAVinv))))
- if(options_.ACES_solver)
- disp('ERROR! saving PI_gensys_data_dump');
- save PI_gensys_data_dump
- error('PI_gensys: Inversion of poss. zero matrix UAVinv=inv(U02''*a1*V02)!');
- else
- warning('PI_gensys: Evading inversion of zero matrix UAVinv=inv(U02''*a1*V02)!');
- eu=[0,0];
- return
- end
+ warning('PI_gensys: Evading inversion of zero matrix UAVinv=inv(U02''*a1*V02)!');
+ eu=[0,0];
+ return
end
-catch
- errmsg=lasterror;
- warning(['error callig PI_gensys_singularC: ' errmsg.message ],'errmsg.identifier');
- %error('errcode',['error callig PI_gensys_singularC: ' errmsg.message ]);
+catch ME
+ warning(['error callig PI_gensys_singularC: ' ME.message ],'errmsg.identifier');
end
%
% Define TT1, TT2
@@ -128,7 +120,6 @@ G21=zeros(FL_RANK,(n-FL_RANK));
G22=zeros(FL_RANK,FL_RANK);
G23=eye(FL_RANK);
%H2=zeros(FL_RANK,NX);
-num_inst=0;
% New Definitions
Ze11=zeros(NX,NX);
@@ -146,34 +137,6 @@ G1pi=[Ze11 Ze12 Ze134 Ze134; P1 G11 G12 G13; Ze31 G21 G22 G23; P3 G31 G32 G33];
impact=[eye(NX,NX); zeros(n+FL_RANK,NX)];
-if(options_.ACES_solver)
- if isfield(lq_instruments,'names')
- num_inst=size(lq_instruments.names,1);
- if num_inst>0
- i_var=lq_instruments.inst_var_indices;
- N1=UAVinv*U02'*lq_instruments.B1;
- N3=-FF*N1+Sinv*U01'*lq_instruments.B1;
- else
- error('WARNING: There are no instrumnets for ACES!');
- end
- lq_instruments.N1=N1;
- lq_instruments.N3=N3;
- else
- error('WARNING: There are no instrumnets for ACES!');
- end
- E3=V02*[P1 G11 G12 G13];
- E3= E3+ [zeros(size(V01,1),size(E3,2)-size(V01,2)) V01];
- E2=-E3;
- E5=-V02*N1;
- DD=[zeros(NX,size(N1,2));N1; zeros(FL_RANK,size(N1,2));N3];
- II=eye(num_inst);
- GAMMA=[ E3 -E5 %zeros(size(E3,1),num_inst);
- zeros(num_inst,size(E3,2)), II;
- ];
- eu =[1; 1], nmat=[], gev=[];
- return % do not check B&K compliancy
-end
-
G0pi=eye(n+FL_RANK+NX);
try
if isoctave && octave_ver_less_than('9')
@@ -182,19 +145,12 @@ try
else
[a, b, q, z]=qz(G0pi,G1pi);
end
-catch
+catch ME
try
- lerror=lasterror;
- disp(['PI_Gensys: ' lerror.message]);
+ disp(['PI_Gensys: ' ME.message]);
if 0==strcmp('MATLAB:qz:matrixWithNaNInf',lerror.identifier)
disp '** Unexpected Error PI_Gensys:qz: ** :';
- button=questdlg('Continue Y/N?','Unexpected Error in qz','No','Yes','Yes');
- switch button
- case 'No'
- error ('Terminated')
- %case 'Yes'
-
- end
+ error ('Unexpected Error in qz: Terminated')
end
G1pi=[];impact=[];nmat=[]; gev=[];
eu=[-2;-2];
@@ -237,7 +193,7 @@ for i=1:nn
zxz=1;
end
end
-div ;
+
if ~zxz
[a, b, ~, z]=qzdiv(div,a,b,q,z);
end
@@ -251,7 +207,7 @@ if zxz
nmat=[]; %;gev=[]
return
end
-if (FL_RANK ~= nunstab && ~options_.ACES_solver)
+if FL_RANK ~= nunstab
disp(['Number of unstable variables ' num2str(nunstab)]);
disp( ['does not match number of expectational equations ' num2str(FL_RANK)]);
nmat=[];% gev=[];
@@ -260,7 +216,6 @@ if (FL_RANK ~= nunstab && ~options_.ACES_solver)
end
% New Definitions
-z1=z(:,1:n+NX)';
z2=z(:,n+NX+1:n+NX+FL_RANK)';
% New N Matrix by J Pearlman
diff --git a/matlab/partial_information/PI_gensys_singularC.m b/matlab/partial_information/PI_gensys_singularC.m
index 15a2a6f75ccee561fcbb43b4f96529a305bb9e8a..73e985a679c14f9bc4f91a68ba11366b2035d95a 100644
--- a/matlab/partial_information/PI_gensys_singularC.m
+++ b/matlab/partial_information/PI_gensys_singularC.m
@@ -28,8 +28,8 @@ level=level+1
if level>100
error( ' PI_gensys_singularC recurssion exceeeded its maximum of 100 iterations! ');
end
-warning('', '');
-M1=[];M2=[]; UAVinv=[];
+
+UAVinv=[];
%
% Find SVD of a0, and create partitions of U, S and V
%
diff --git a/matlab/partial_information/add_auxiliary_variables_to_steadystate.m b/matlab/partial_information/add_auxiliary_variables_to_steadystate.m
deleted file mode 100644
index 257d055bb255c6d89fd271a76d35bd062af35616..0000000000000000000000000000000000000000
--- a/matlab/partial_information/add_auxiliary_variables_to_steadystate.m
+++ /dev/null
@@ -1,41 +0,0 @@
-function ys1 = add_auxiliary_variables_to_steadystate(ys,aux_vars,fname, ...
- exo_steady_state, exo_det_steady_state,params, byte_code)
-% Add auxiliary variables to the steady state vector
-
-% Copyright © 2009-2024 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-
-global M_ options_
-
-n = length(aux_vars);
-ys1 = [ys;zeros(n,1)];
-
-for i=1:n+1
- if byte_code
- res = bytecode('static','evaluate', M_, options_, ys1,...
- [exo_steady_state; ...
- exo_det_steady_state],params);
- else
- res = feval([fname '.sparse.static_resid'], ys1,...
- [exo_steady_state; ...
- exo_det_steady_state],params);
- end
- for j=1:n
- el = aux_vars(j).endo_index;
- ys1(el) = ys1(el)-res(el);
- end
-end
diff --git a/matlab/partial_information/disc_riccati_fast.m b/matlab/partial_information/disc_riccati_fast.m
index dabc2a6531129ba4f643579dd3f18e0c1c44cc0f..2ca75fdc23f39ca5e125a140d001721aab2fb355 100644
--- a/matlab/partial_information/disc_riccati_fast.m
+++ b/matlab/partial_information/disc_riccati_fast.m
@@ -11,7 +11,7 @@ function Z=disc_riccati_fast(F,D,R,H,ch)
% V.2 22/10/06
% =================================================================
-% Copyright © 2006-2017 Dynare Team
+% Copyright © 2006-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -34,7 +34,6 @@ else
flag_ch = 1;
end
-
% intialisation
tol = 1e-10; % iteration convergence threshold
P0=H;
@@ -43,7 +42,7 @@ if ~any(R) % i.e. ==0
warning('Dangerously evading inversion of zero matrix!');
Y0=zeros(size(R));
else
- Y0=D'*inv(R)*D;
+ Y0=D'/R*D;
end
POYO=P0*Y0;
I=eye(size(POYO));
@@ -63,7 +62,6 @@ while matd > tol && count < 100
X0=X1;
Y0=Y1;
count=count+1;
- % matd;
end
Z=(P0+P0')/2;
@@ -72,15 +70,7 @@ Z=(P0+P0')/2;
if count==100
matd
error('Riccati not converged fast enough!');
- % error.identifier='Riccati not converged!'
- % error
end
-%if count >5
-% disp('Riccati count= ');
-% count
-%end
-
-clear X0 X1 Y0 Y1 P1 I INVPY;
% Check that X is positive definite
if flag_ch==1
diff --git a/matlab/partial_information/dr1_PI.m b/matlab/partial_information/dr1_PI.m
index b583d61d98190a32f9fd16d66cedd2da64e6ccd9..2a3b289fe57c82d0f30325d0f50dd7e75470f13d 100644
--- a/matlab/partial_information/dr1_PI.m
+++ b/matlab/partial_information/dr1_PI.m
@@ -1,5 +1,5 @@
-function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
-% function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
+function [dr,info,oo_] = dr1_PI(dr,M_,options_,oo_)
+% function [dr,info,oo_] = dr1_PI(dr,M_,options_,oo_)
% Computes the reduced form solution of a rational expectation model first
% order
% approximation of the Partial Information stochastic model solver around the deterministic steady state).
@@ -14,8 +14,6 @@ function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
%
% INPUTS
% dr [matlab structure] Decision rules for stochastic simulations.
-% task [integer] if task = 0 then dr1 computes decision rules.
-% if task = 1 then dr1 computes eigenvalues.
% M_ [matlab structure] Definition of the model.
% options_ [matlab structure] Global options.
% oo_ [matlab structure] Results
@@ -29,8 +27,6 @@ function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
% info=4: BK order condition not satisfied info(2) contains "distance"
% indeterminacy.
% info=5: BK rank condition not satisfied.
-% M_ [matlab structure]
-% options_ [matlab structure]
% oo_ [matlab structure]
%
% ALGORITHM
@@ -57,7 +53,6 @@ function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global lq_instruments;
info = 0;
options_ = set_default_option(options_,'qz_criterium',1.000001);
@@ -65,166 +60,31 @@ options_ = set_default_option(options_,'qz_criterium',1.000001);
xlen = M_.maximum_endo_lead + M_.maximum_endo_lag + 1;
if options_.aim_solver
- options_.aim_solver = false;
- warning('You can not use AIM with Part Info solver. AIM ignored');
-end
+ error('Anderson and Moore AIM solver is not compatible with Partial Information models');
+end % end if useAIM and...
+
if (options_.order > 1)
warning('You can not use order higher than 1 with Part Info solver. Order 1 assumed');
options_.order =1;
end
-% expanding system for Optimal Linear Regulator
-if options_.ramsey_policy && ~options_.ACES_solver
- if isfield(M_,'orig_model')
- orig_model = M_.orig_model;
- M_.endo_nbr = orig_model.endo_nbr;
- M_.endo_names = orig_model.endo_names;
- M_.lead_lag_incidence = orig_model.lead_lag_incidence;
- M_.maximum_lead = orig_model.maximum_lead;
- M_.maximum_endo_lead = orig_model.maximum_endo_lead;
- M_.maximum_lag = orig_model.maximum_lag;
- M_.maximum_endo_lag = orig_model.maximum_endo_lag;
- end
- o_jacobian_flag = options_.jacobian_flag;
- options_.jacobian_flag = false;
- oo_.steady_state = dynare_solve('ramsey_static', oo_.steady_state, ...
- options_.ramsey.maxit, options_.solve_tolf, options_.solve_tolx, ...
- options_, M_, options_, oo_, it_);
- options_.jacobian_flag = o_jacobian_flag;
- [~,~,multbar] = ramsey_static(oo_.steady_state,M_,options_,oo_,it_);
- [jacobia_,M_] = ramsey_dynamic(oo_.steady_state,multbar,M_,options_,oo_,it_);
- klen = M_.maximum_lag + M_.maximum_lead + 1;
- dr.ys = [oo_.steady_state;zeros(M_.exo_nbr,1);multbar];
-else
- klen = M_.maximum_lag + M_.maximum_lead + 1;
- iyv = M_.lead_lag_incidence';
- iyv = iyv(:);
- iyr0 = find(iyv) ;
- it_ = M_.maximum_lag + 1 ;
-
- if M_.exo_nbr == 0
- oo_.exo_steady_state = [] ;
- end
-
+klen = M_.maximum_lag + M_.maximum_lead + 1;
+iyv = M_.lead_lag_incidence';
- if options_.ACES_solver
- sim_ruleids=[];
- tct_ruleids=[];
- if size(M_.equations_tags,1)>0 % there are tagged equations, check if they are aceslq rules
- for teq=1:size(M_.equations_tags,1)
- if strcmp(M_.equations_tags(teq,3),'aceslq_sim_rule')
- sim_ruleids=[sim_ruleids cell2mat(M_.equations_tags(teq,1))]
- end
- if strcmp(M_.equations_tags(teq,3),'aceslq_tct_rule')
- tct_ruleids=[tct_ruleids cell2mat(M_.equations_tags(teq,1))]
- end
- end
- end
- lq_instruments.sim_ruleids=sim_ruleids;
- lq_instruments.tct_ruleids=tct_ruleids;
- %if isfield(lq_instruments,'xsopt_SS') %% changed by BY
- [~, lq_instruments.xsopt_SS,lq_instruments.lmopt_SS,s2,check] = opt_steady_get;%% changed by BY
- [~, DYN_Q] = QPsolve(lq_instruments, s2, check); %% added by BY
- z = repmat(lq_instruments.xsopt_SS,1,klen);
- else
- z = repmat(dr.ys,1,klen);
- end
- z = z(iyr0) ;
- [~,jacobia_] = feval([M_.fname '.dynamic'],z,[oo_.exo_simul ...
- oo_.exo_det_simul], M_.params, dr.ys, it_);
+if M_.exo_nbr == 0
+ oo_.exo_steady_state = [] ;
+end
- if options_.ACES_solver && (length(sim_ruleids)>0 || length(tct_ruleids)>0 )
- if length(sim_ruleids)>0
- sim_rule=jacobia_(sim_ruleids,:);
- % uses the subdirectory - BY
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_sim_rule.txt'], 'sim_rule', '-ascii', '-double', '-tabs');
- end
- if length(tct_ruleids)>0
- tct_rule=jacobia_(tct_ruleids,:);
- % uses the subdirectory - BY
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_tct_rule.txt'], 'tct_rule', '-ascii', '-double', '-tabs');
- end
- aces_ruleids=union(tct_ruleids,sim_ruleids);
- j_size=size(jacobia_,1);
- j_rows=1:j_size;
- j_rows = setxor(j_rows,aces_ruleids);
- jacobia_=jacobia_(j_rows ,:);
- end
-end
+z = repmat(dr.ys,1,klen);
+z = z(find(iyv(:))) ;
+[~,jacobia_] = feval([M_.fname '.dynamic'],z,[oo_.exo_simul ...
+ oo_.exo_det_simul], M_.params, dr.ys, M_.maximum_lag + 1);
if options_.debug
save([M_.dname filesep 'Output' filesep M_.fname '_debug.mat'],'jacobia_')
end
-dr=set_state_space(dr,M_);
-nstatic = M_.nstatic;
-nfwrd = M_.nfwrd;
-nspred = M_.nspred;
-nboth = M_.nboth;
-order_var = dr.order_var;
-nd = M_.nspred+M_.nsfwrd;
-nz = nnz(M_.lead_lag_incidence);
-
-sdyn = M_.endo_nbr - nstatic;
-
-k0 = M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var);
-k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
-
-if options_.aim_solver
- error('Anderson and Moore AIM solver is not compatible with Partial Information models');
-end % end if useAIM and...
-
-% If required, try PCL86 solver, that is, if not the check being
-% performed only and if it is 1st order
-% create sparse, extended jacobia AA:
-nendo=M_.endo_nbr; % = size(aa,1)
-
-
-if(options_.ACES_solver)
- %if ~isfield(lq_instruments,'names')
- if isfield(options_,'instruments')
- lq_instruments.names=options_.instruments;
- end
- %end
- if isfield(lq_instruments,'names')
- num_inst=size(lq_instruments.names,1);
- if ~isfield(lq_instruments,'inst_var_indices') && num_inst>0
- for i=1:num_inst
- i_tmp = strmatch(deblank(lq_instruments.names(i,:)), M_.endo_names,'exact');
- if isempty(i_tmp)
- error ('One of the specified instrument variables does not exist') ;
- else
- i_var(i) = i_tmp;
- end
- end
- lq_instruments.inst_var_indices=i_var;
- elseif size(lq_instruments.inst_var_indices)>0
- i_var=lq_instruments.inst_var_indices;
- if ~num_inst
- num_inst=size(lq_instruments.inst_var_indices);
- end
- else
- i_var=[];
- num_inst=0;
- end
- if size(i_var,2)>0 && size(i_var,2)==num_inst
- m_var=zeros(nendo,1);
- for i=1:nendo
- if isempty(find(i_var==i))
- m_var(i)=i;
- end
- end
- m_var=nonzeros(m_var);
- lq_instruments.m_var=m_var;
- else
- error('WARNING: There are no instrumnets for ACES!');
- end
- else %if(options_.ACES_solver==1)
- error('WARNING: There are no instrumnets for ACES!');
- end
-end
-
% find size xlen of the state vector Y and of A0, A1 and A2 transition matrices:
% it is the sum the all i variables's lag/lead representations,
% for each variable i representation being defined as:
@@ -242,75 +102,36 @@ end
% partition jacobian:
jlen=M_.nspred+M_.nsfwrd+M_.endo_nbr+M_.exo_nbr; % length of jacobian
-PSI=-jacobia_(:, jlen-M_.exo_nbr+1:end); % exog
- % first transpose M_.lead_lag_incidence';
+% first transpose M_.lead_lag_incidence';
lead_lag=M_.lead_lag_incidence';
-max_lead_lag=zeros(nendo,2); % lead/lag representation in Y for each endogenous variable i
if ( M_.maximum_lag <= 1) && (M_.maximum_lead <= 1)
xlen=size(jacobia_,1);%nendo;
AA0=zeros(xlen,xlen); % empty A0
AA2=AA0; % empty A2 and A3
AA3=AA0;
- if xlen==nendo % && M_.maximum_lag <=1 && M_.maximum_lead <=1 % apply a shortcut
- AA1=jacobia_(:,nspred+1:nspred+nendo);
+ if xlen==M_.endo_nbr % && M_.maximum_lag <=1 && M_.maximum_lead <=1 % apply a shortcut
+ AA1=jacobia_(:,M_.nspred+1:M_.nspred+M_.endo_nbr);
if M_.maximum_lead ==1
fnd = find(lead_lag(:,M_.maximum_lag+2));
AA0(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,M_.maximum_lag+2))); %forwd jacobian
end
- if nspred>0 && M_.maximum_lag ==1
+ if M_.nspred>0 && M_.maximum_lag ==1
fnd = find(lead_lag(:,1));
AA2(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,1))); %backward
end
- elseif options_.ACES_solver % more endo vars than equations in jacobia_
- if nendo-xlen==num_inst
- PSI=[PSI;zeros(num_inst, M_.exo_nbr)];
- % AA1 contemporary
- AA_all=jacobia_(:,nspred+1:nspred+nendo);
- AA1=AA_all(:,lq_instruments.m_var); % endo without instruments
- lq_instruments.ij1=AA_all(:,lq_instruments.inst_var_indices); % instruments only
- lq_instruments.B1=-[lq_instruments.ij1; eye(num_inst)];
- AA1=[AA1, zeros(xlen,num_inst); zeros(num_inst,xlen), eye(num_inst)];
- %PSI=[PSI; zeros(num_inst,M_.exo_nbr)];
- if M_.maximum_lead ==1 % AA0 forward looking
- AA_all(:,:)=0.0;
- fnd = find(lead_lag(:,M_.maximum_lag+2));
- AA_all(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,M_.maximum_lag+2))); %forwd jacobian
- AA0=AA_all(:,lq_instruments.m_var);
- lq_instruments.ij0=AA_all(:,lq_instruments.inst_var_indices); % instruments only
- lq_instruments.B0=[lq_instruments.ij0; eye(num_inst)];
- AA0=[AA0, zeros(xlen,num_inst); zeros(num_inst,xlen+num_inst)];
- end
- if nspred>0 && M_.maximum_lag ==1
- AA_all(:,:)=0.0;
- fnd = find(lead_lag(:,1));
- AA_all(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,1))); %backward
- AA2=AA_all(:,lq_instruments.m_var);
- lq_instruments.ij2=AA_all(:,lq_instruments.inst_var_indices); % instruments only
- lq_instruments.B2=[lq_instruments.ij2; eye(num_inst)];
- AA2=[AA2, lq_instruments.ij2 ; zeros(num_inst,xlen+num_inst)];
- end
- else
- error('ACES number of instruments does match');
- end
else
error('More than one lead or lag in the jabian');
end
- if M_.orig_endo_nbr<nendo
- % findif there are any expecatations at time t
+ if M_.orig_endo_nbr<M_.endo_nbr
+ % findif there are any expectations at time t
exp_0= strmatch('AUX_EXPECT_LEAD_0_', M_.endo_names);
- num_exp_0=size(exp_0,1);
- if num_exp_0>0
- AA3(:,exp_0)=AA1(:,exp_0);
- XX0=zeros(nendo,num_exp_0);
- AA1(:,exp_0)=XX0(:,1:num_exp_0)
+ if ~isempty(exp_0)
+ error('Partial Information does not support EXPECTATION(0) operators')
end
end
end
-PSI=-[[zeros(xlen-nendo,M_.exo_nbr)];[jacobia_(:, jlen-M_.exo_nbr+1:end)]]; % exog
+PSI=-[[zeros(xlen-M_.endo_nbr,M_.exo_nbr)];[jacobia_(:, jlen-M_.exo_nbr+1:end)]]; % exog
cc=0;
-NX=M_.exo_nbr; % no of exogenous varexo shock variables.
-NETA=nfwrd+nboth; % total no of exp. errors set to no of forward looking equations
-FL_RANK=rank(AA0); % nfwrd+nboth; % min total no of forward looking equations and vars
try
% call [G1pi,C,impact,nmat,TT1,TT2,gev,eu]=PI_gensys(a0,a1,a2,c,PSI,NX,NETA,NO_FL_EQS)
@@ -323,32 +144,15 @@ try
% (b) matrices TT1, TT2 that relate y(t) to these states:
% y(t)=[TT1 TT2][s(t)' x(t)']'.
- if(options_.ACES_solver)
- if isfield(lq_instruments,'xsopt_SS')
- SSbar= diag([lq_instruments.xsopt_SS(m_var)]);% lq_instruments.xsopt_SS(lq_instruments.inst_var_indices)]);
- insSSbar=repmat(lq_instruments.xsopt_SS(lq_instruments.inst_var_indices)',nendo-num_inst,1);
- else
- SSbar= diag([dr.ys(m_var)]);%; dr.ys(lq_instruments.inst_var_indices)]);%(oo_.steady_state);
- insSSbar=repmat(dr.ys(lq_instruments.inst_var_indices)',nendo-num_inst,1);
- end
- SSbar=diag([diag(SSbar);diag(eye(num_inst))]);
- insSSbar=[insSSbar;diag(eye(num_inst))];
-
- AA0=AA0*SSbar;
- AA1=AA1*SSbar;
- AA2=AA2*SSbar;
- lq_instruments.B1=(lq_instruments.B1).*insSSbar;
- end
%% for expectational models when complete
if any(AA3)
- AA3=AA3*SSbar;
- [G1pi,CC,impact,nmat,TT1,TT2,gev,eu, DD, E2,E5, GAMMA, FL_RANK]=PI_gensysEXP(AA0,AA1,-AA2,AA3,cc,PSI,NX,NETA,FL_RANK, M_, options_);
+ error('Unsupported case')
else
- [G1pi,CC,impact,nmat,TT1,TT2,gev,eu, DD, E2,E5, GAMMA, FL_RANK]=PI_gensys(AA0,AA1,-AA2,AA3,cc,PSI,NX,NETA,FL_RANK, M_, options_);
+ [G1pi,CC,impact,nmat,TT1,TT2,gev,eu, DD, E2,E5, GAMMA, FL_RANK]=PI_gensys(AA0,AA1,-AA2,cc,PSI,M_.exo_nbr);
end
% reuse some of the bypassed code and tests that may be needed
- if (eu(1) ~= 1 || eu(2) ~= 1) && ~options_.ACES_solver
+ if (eu(1) ~= 1 || eu(2) ~= 1)
info(1) = abs(eu(1)+eu(2));
info(2) = 1.0e+8;
% return
@@ -363,91 +167,12 @@ try
dr.PI_gev=gev;
dr.PI_eu=eu;
dr.PI_FL_RANK=FL_RANK;
- %dr.ys=zeros(nendo); % zero steady state
dr.ghx=G1pi;
dr.ghu=impact;
dr.eigval = eig(G1pi);
dr.rank=FL_RANK;
-
- if options_.ACES_solver
- betap=options_.planner_discount;
- sigma_cov=M_.Sigma_e;
- % get W - BY
- W=(1-betap)*GAMMA'*DYN_Q*GAMMA;
- %W=[0]
- ACES.A=G1pi;
- ACES.C=impact; % (:,1);
- ACES.D=DD; %=impact (:,20);
- ACES.E2=E2;
- ACES.E5=E5;
- ACES.GAMMA=GAMMA;
- ACES_M=size(G1pi,2)-FL_RANK;
- ACES_NM=FL_RANK;
- ACES.M=ACES_M;
- ACES.NM=FL_RANK;
- % added by BY
- ACES.Q=DYN_Q;
- ACES.W=W;
- NY=nendo-num_inst;
-
- % save the followings in a subdirectory - BY
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_Matrices'], 'ACES');
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_GAMMA'], 'GAMMA');
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_A.txt'], 'G1pi', '-ascii', '-double', '-tabs');
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_C.txt'], 'impact','-ascii', '-double', '-tabs');
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_D.txt'], 'DD', '-ascii', '-double', '-tabs');
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_E2.txt'], 'E2', '-ascii', '-double', '-tabs');
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_E5.txt'], 'E5', '-ascii', '-double', '-tabs');
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_GAMMA.txt'], 'GAMMA', '-ascii', '-double', '-tabs');
- %save ([M_.fname '_ACESLQ_M.txt'], 'ACES_M', '-ascii', '-tabs');
- %save ([M_.fname '_ACESLQ_NM.txt'], 'ACES_NM', '-ascii', '-tabs');
- %save ([M_.fname '_ACESLQ_betap.txt'], 'betap', '-ascii', '-tabs');
- %save ([M_.fname '_ACESLQ_NI.txt'], 'num_inst', '-ascii', '-tabs');
- %save ([M_.fname '_ACESLQ_ND.txt'], 'NX', '-ascii', '-tabs');
- %save ([M_.fname '_ACESLQ_NY.txt'], 'NY', '-ascii', '-tabs');
- ACES_VARS=char(M_.endo_names);
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_VARS.txt'], 'ACES_VARS', '-ascii', '-tabs');
- % added by BY
- % save the char array ACES_VARS into .txt as it is
- fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_VARnames.txt'),'wt');
- ACES_VARS =[ACES_VARS repmat(sprintf('\n'),size(ACES_VARS,1),1)];
- fwrite(fid,ACES_VARS.');
- fclose(fid);
- % save as integers
- fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_M.txt'),'wt');
- fprintf(fid,'%d\n',ACES_M);
- fclose(fid);
- fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_NM.txt'),'wt');
- fprintf(fid,'%d\n',ACES_NM);
- fclose(fid);
- fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_betap.txt'),'wt');
- fprintf(fid,'%d\n',betap);
- fclose(fid);
- fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_NI.txt'),'wt');
- fprintf(fid,'%d\n',num_inst);
- fclose(fid);
- fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_ND.txt'),'wt');
- fprintf(fid,'%d\n',NX);
- fclose(fid);
- fid = fopen(strcat(ACES_DirectoryName,'/',M_.fname,'_ACESLQ_NY.txt'),'wt');
- fprintf(fid,'%d\n',NY);
- fclose(fid);
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_Q.txt'], 'DYN_Q', '-ascii', '-double', '-tabs');
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_W.txt'], 'W', '-ascii', '-double', '-tabs');
- save ([ACES_DirectoryName,'/',M_.fname '_ACESLQ_SIGMAE.txt'], 'sigma_cov', '-ascii', '-double', '-tabs');
- end
-
-catch
- lerror=lasterror;
- if options_.ACES_solver
- disp('Problem with using Part Info ACES solver:');
- error(lerror.message);
- else
- disp('Problem with using Part Info solver');
- error(lerror.message);
- end
-end
-
-% TODO:
-% if options_.loglinear == 1
-% if exogenous deterministic variables
+
+catch ME
+ disp('Problem with using Part Info solver');
+ rethrow(ME);
+end
\ No newline at end of file
diff --git a/matlab/stochastic_solver/stoch_simul.m b/matlab/stochastic_solver/stoch_simul.m
index 3d0cc4e0fc4c2a2b4b2205c5a5b1341a8c191c1a..376de52368289594b0e0a1abfa3f4275cd53a234 100644
--- a/matlab/stochastic_solver/stoch_simul.m
+++ b/matlab/stochastic_solver/stoch_simul.m
@@ -63,10 +63,10 @@ if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;
end
-if options_.partial_information || options_.ACES_solver
+if options_.partial_information
PI_PCL_solver = 1;
if options_.order ~= 1
- warning('stoch_simul:: forcing order=1 since you are using partial_information or ACES solver')
+ warning('stoch_simul:: forcing order=1 since you are using partial_information')
options_.order = 1;
end
else
@@ -93,7 +93,7 @@ check_model(M_);
oo_.dr=set_state_space(oo_.dr,M_);
if PI_PCL_solver
- [oo_.dr, info] = PCL_resol(oo_.steady_state,0);
+ [oo_.dr, info] = PCL_resol(M_, options_, oo_);
elseif options_.discretionary_policy
if ~options_.order==1
error('discretionary_policy: only linear-quadratic problems can be solved');
@@ -201,7 +201,7 @@ end
if ~options_.nomoments
if PI_PCL_solver
- PCL_Part_info_moments(0, PCL_varobs, oo_.dr, i_var);
+ oo_=PCL_Part_info_moments(M_, oo_, options_, PCL_varobs, oo_.dr, i_var);
elseif options_.periods == 0
if options_.order == 1 || (options_.order == 2 && ~options_.pruning)
oo_=disp_th_moments(oo_.dr,var_list,M_,options_,oo_);
diff --git a/meson.build b/meson.build
index 661ef4731ac0e53f854243b4f5d74a60d1c70248..a8d08c729704f2558fe75d14263e164fb1a3b4a3 100644
--- a/meson.build
+++ b/meson.build
@@ -1129,6 +1129,8 @@ mod_and_m_tests = [
{ 'test' : [ 'stochastic_simulations/fs2000_k_order_simul.mod' ] },
{ 'test' : [ 'partial_information/PItest3aHc0PCLsimModPiYrVarobsAll.mod' ] },
{ 'test' : [ 'partial_information/PItest3aHc0PCLsimModPiYrVarobsCNR.mod' ] },
+ { 'test' : [ 'partial_information/Partial_ramsey.mod' ] },
+ { 'test' : [ 'partial_information/Partial_steady_state_model.mod' ] },
{ 'test' : [ 'arima/mod1.mod',
'arima/mod1a.mod',
'arima/mod1b.mod',
diff --git a/tests/partial_information/Partial_ramsey.mod b/tests/partial_information/Partial_ramsey.mod
new file mode 100644
index 0000000000000000000000000000000000000000..6912c172f32b99754607e9657709465964084a9a
--- /dev/null
+++ b/tests/partial_information/Partial_ramsey.mod
@@ -0,0 +1,177 @@
+/*
+ * This file implements the optimal monetary policy under commitment exercise
+ * in Jordi Galà (2008): Monetary Policy, Inflation, and the Business Cycle,
+ * Princeton University Press, Chapter 5.1.2
+ *
+ * It demonstrates how to use the ramsey_model command of Dynare.
+ *
+ * Notes:
+ * - all model variables are expressed in deviations from steady state, i.e.
+ * in contrast to to the chapter, both the nominal interest rate and
+ * natural output are not in log-levels, but rather mean 0
+ *
+ * This implementation was written by Johannes Pfeifer. In case you spot mistakes,
+ * email me at jpfeifer@gmx.de
+ *
+ * Please note that the following copyright notice only applies to this Dynare
+ * implementation of the model.
+ */
+
+/*
+ * Copyright 2015-2024 Johannes Pfeifer
+ *
+ * This is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * It is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * For a copy of the GNU General Public License,
+ * see <https://www.gnu.org/licenses/>.
+ */
+
+
+var pi ${\pi}$ (long_name='inflation')
+ y_gap ${tilde y}$ (long_name='output gap')
+ y_nat ${y^{nat}}$ (long_name='natural output')
+ y ${y}$ (long_name='output')
+ r_e ${r^{e}}$ (long_name='efficient interest rate')
+ y_e ${y^{nat}}$ (long_name='efficient output')
+ x ${x}$ (long_name='welfare-relevant output gap')
+ r_nat ${r^{nat}}$ (long_name='natural interest rate')
+ r_real ${r^r}$ (long_name='real interest rate')
+ i ${i}$ (long_name='nominal interest rate')
+ n ${n}$ (long_name='hours worked')
+ m_growth_ann ${\Delta m}$ (long_name='money growth')
+ u ${u}$ (long_name='AR(1) cost push shock process')
+ a ${a}$ (long_name='AR(1) technology shock process')
+ r_real_ann ${r^{r,ann}}$ (long_name='annualized real interest rate')
+ i_ann ${i^{ann}}$ (long_name='annualized nominal interest rate')
+ r_nat_ann ${r^{nat,ann}}$ (long_name='annualized natural interest rate')
+ pi_ann ${\pi^{ann}}$ (long_name='annualized inflation rate')
+ p ${p}$ (long_name='price level')
+ ;
+
+varexo eps_a ${\varepsilon_a}$ (long_name='technology shock')
+ eps_u ${\varepsilon_u}$ (long_name='monetary policy shock');
+
+parameters alppha ${\alppha}$ (long_name='capital share')
+ betta ${\beta}$ (long_name='discount factor')
+ rho_a ${\rho_a}$ (long_name='autocorrelation technology shock')
+ rho_u ${\rho_{u}}$ (long_name='autocorrelation cost push shock')
+ siggma ${\sigma}$ (long_name='log utility')
+ phi ${\phi}$ (long_name='unitary Frisch elasticity')
+ phi_y ${\phi_{y}}$ (long_name='output feedback Taylor Rule')
+ eta ${\eta}$ (long_name='semi-elasticity of money demand')
+ epsilon ${\epsilon}$ (long_name='demand elasticity')
+ theta ${\theta}$ (long_name='Calvo parameter')
+ ;
+%----------------------------------------------------------------
+% Parametrization, p. 52
+%----------------------------------------------------------------
+siggma = 1;
+phi=1;
+phi_y = .5/4;
+theta=2/3;
+rho_u = 0;
+rho_a = 0.9;
+betta = 0.99;
+eta =4;
+alppha=1/3;
+epsilon=6;
+
+
+
+%----------------------------------------------------------------
+% First Order Conditions
+%----------------------------------------------------------------
+
+model(linear);
+//Composite parameters
+#Omega=(1-alppha)/(1-alppha+alppha*epsilon); //defined on page 47
+#psi_n_ya=(1+phi)/(siggma*(1-alppha)+phi+alppha); //defined on page 48
+#lambda=(1-theta)*(1-betta*theta)/theta*Omega; //defined on page 47
+#kappa=lambda*(siggma+(phi+alppha)/(1-alppha)); //defined on page 49
+#alpha_x=kappa/epsilon; //defined on page 96
+#phi_pi=(1-rho_u)*kappa*siggma/(alpha_x)+rho_u; //defined on page 101
+
+//1. Definition efficient interest rate, below equation (4)
+r_e=siggma*(y_e(+1)-y_e);
+
+//2. Definition efficient output
+y_e=psi_n_ya*a;
+
+//3. Definition linking various output gaps, bottom page 96
+y_gap=x+(y_e-y_nat);
+
+//4. New Keynesian Phillips Curve eq. (2)
+pi=betta*pi(+1)+kappa*x + u;
+
+//5. Dynamic IS Curve eq. (4)
+x=-1/siggma*(i-pi(+1)-r_e)+x(+1);
+
+//6. Definition natural rate of interest eq. (23)
+r_nat=siggma*psi_n_ya*(a(+1)-a);
+
+//7. Definition real interest rate
+r_real=i-pi(+1);
+
+//8. Natural output
+y_nat=phi*a;
+
+//9. Definition output gap
+y_gap=y-y_nat;
+
+//10. cost push shock, equation (3)
+u=rho_u*u(-1)+eps_u;
+
+//11. TFP shock
+a=rho_a*a(-1)+eps_a;
+
+//12. Production function (eq. 13)
+y=a+(1-alppha)*n;
+
+//13. Money growth (derived from eq. (4))
+m_growth_ann=4*(y-y(-1)-eta*(i-i(-1))+pi);
+
+//14. Annualized nominal interest rate
+i_ann=4*i;
+
+//15. Annualized real interest rate
+r_real_ann=4*r_real;
+
+//16. Annualized natural interest rate
+r_nat_ann=4*r_nat;
+
+//17. Annualized inflation
+pi_ann=4*pi;
+
+//18. Definition price level
+pi=p-p(-1);
+
+//19. Interest Rate Rule that implements optimal solution, eq. (10)
+% i=r_e+phi_pi*pi;
+end;
+
+%----------------------------------------------------------------
+% define shock variances
+%---------------------------------------------------------------
+
+shocks;
+var eps_u = 1;
+var eps_a = 1;
+end;
+
+//planner objective using alpha_x expressed as function of deep parameters
+planner_objective pi^2 +(((1-theta)*(1-betta*theta)/theta*((1-alppha)/(1-alppha+alppha*epsilon)))*(siggma+(phi+alppha)/(1-alppha)))/epsilon*y_gap^2;
+
+ramsey_model(instruments=(i),planner_discount=betta);
+
+varobs pi x y u a;
+stoch_simul(order=1,irf=13,partial_information) x pi u;
+
+check;
\ No newline at end of file
diff --git a/tests/partial_information/Partial_steady_state_model.mod b/tests/partial_information/Partial_steady_state_model.mod
new file mode 100644
index 0000000000000000000000000000000000000000..cd4d6d9bc214a815ef847c11b92b00ec619acf0a
--- /dev/null
+++ b/tests/partial_information/Partial_steady_state_model.mod
@@ -0,0 +1,183 @@
+var r, pi, n, w, y, mc, x, rr, yobs, piobs,robs;
+
+var y_star, n_star, w_star;
+var y_f,pi_f,y_e,pi_e, pi_exp,y_exp;
+%Exogenous process (common to both specifications)
+var a, nu, u;
+%Shocks (common to both specifications)
+varexo eps_a,eps_nu,eps_u;
+%model parameters (common to both specifications)
+parameters h, gama, phi, bet, var_omega, kappa, iota_r, iota_y, iota_pi,iota_piexp;
+%shock persistence (common to both specifications)
+parameters rho_a,rho_nu,rho_u;
+%some steady state variables (common to both specifications)
+parameters a_ss, n_ss, mc_ss;
+% learning parameters
+parameters theta, rho,sigmaScale;
+
+parameters alpha_pi_e, alpha_pi_f, alpha_y_f, alpha_y_e;
+
+%----------------------------------------------------------------
+% Parametrization
+%----------------------------------------------------------------
+%model parameters
+h=0.5; %%%%%%%%%%to be matched with moments habit (0,1)
+gama=1; %%%%%%%%%%to be matched with moments inverse of the inter. elasticity of subst.
+phi=1;%%%%%%%%%%to be matched with moments %inverse of the Frisch
+bet=0.99; %%%%%%%%%%set
+var_omega=0.75; %%%%%%%%%%to be matched with moments Calvo's probability
+kappa=0.5;%%%%%%%%%%to be matched with moments Past indexation (0,1)
+%labor
+n_ss=0.2;
+%elasticity of substituion across varieties sigma
+sigma=6;
+%marginal cost
+mc_ss=(sigma-1)/sigma;
+%profit_share
+ profit_share=1/sigma;
+ wage_share=1-profit_share;
+%use the demand for labor
+a_ss=(wage_share)/(mc_ss*n_ss);
+% % Taylor rule parameters
+iota_r=0.2; %set %%%%%%%%%% to be estimated
+iota_y=0.5/4; %set%%%%%%%%%% to be estimated
+iota_pi=1.5; %set%%%%%%%%%% to be estimated
+iota_piexp=0; %set%%%%%%%%%% to be estimated
+rho_a=0.8; %%%%%%%%%%to be matched with moments persistence tech shock
+rho_nu=0.8; %%%%%%%%%%to be matched with moments persistence demand shock
+rho_u=0.5; %%%%%%%%%%to be matched with moments existence monetary policy shock (usually lower than the other two)
+%Behavioral Parameters
+theta=2; %%%%%%%%%% to be estimated
+rho=0.5; %%%%%%%%%% to be estimated
+sigmaScale=0.005; %%%%%%%%%% Scale Parameter of the logistic function
+alpha_pi_e=0.5;
+alpha_pi_f=0.5;
+alpha_y_f=0.5;
+alpha_y_e=0.5;
+
+%----------------------------------------------------------------
+% The Model
+%----------------------------------------------------------------
+model;
+//Composite parameters
+#messy1=((1-var_omega*bet)/(1+bet*kappa))*(1-var_omega)/var_omega; % I coeff. NKPC
+#messy2=bet/(1+bet*kappa); % II coeff. NKPC
+#messy3=kappa/(1+bet*kappa); % III coeff. NKPC
+#epsilon_n=(a_ss*n_ss); %elasticity wrt n
+%sticky prices
+[name='Rep Euler Equation']
+y*(1+h)=y(+1)+h*y(-1)-((1-h)/gama)*(r-pi(+1))+nu;
+%eq2
+[name='Labor supply']
+n=(w-gama*y/(1-h)+gama*h*y(-1)/(1-h))/phi;
+%eq3
+[name='Production function']
+y=epsilon_n*(a+n);
+%eq4
+[name='Labor Demand']
+w=mc+(y-n)+a;
+[name='Rep NKPC']
+pi=messy1*mc+messy2*pi(+1)+messy3*pi(-1);
+%eq9
+[name='Taylor Rule']
+r=(1-iota_r)*(iota_pi*pi+iota_y*(y-y_star)+iota_piexp*pi_exp)+iota_r*r(-1)+u;
+%eq10
+[name='Output Gap']
+x=y-y_star;
+%%%%Flex prices
+%eq11
+[name='Labor Supply - Flex Prices']
+n_star=(w_star-gama*y_star/(1-h)+gama*h*y_star(-1)/(1-h))/phi;
+%eq12
+[name='Production Function - Flex Prices']
+y_star=epsilon_n*(a+n_star);
+%eq13
+[name='Labor Demand - Flex Prices']
+w_star=y_star-n_star+a;
+% % % % % % % % Equation swap
+%eq17
+ [name='Rational Expectation Y']
+ y_f=y(+1);
+% %eq18
+ [name='Rational Expectation PI']
+ pi_f=pi(+1);
+% % Fundamentalist
+%eq19
+[name='Extrapolative Expectation Y']
+y_e=y(-1);
+%eq20
+[name='Extrapolative Expectation PI']
+pi_e=pi(-1);
+%eq21
+[name='Mkt Expectation Y']
+y_exp=alpha_y_f*y_f+alpha_y_e*y_e;
+%eq22
+[name='Mkt Expectation PI']
+pi_exp=alpha_pi_f*pi_f+alpha_pi_e*pi_e;
+%eq31
+[name='Real Interest Rate']
+rr=r-pi(+1);
+% Observation Equations
+yobs=y;
+piobs=pi;
+robs=r;
+%Exogenous process
+%eq33-36
+[name='TFP Shock']
+a=rho_a*a(-1)+eps_a;
+[name='Demand Schock']
+nu=rho_nu*nu(-1)+eps_nu;
+[name='Monetary Policy Shock']
+u=rho_u*u(-1)+eps_u;
+end;
+
+steady_state_model;
+mc=0;
+n=0;
+w=0;
+pi=0;
+piobs=0;
+r=0;
+robs=0;
+y=0;
+yobs=0;
+x=0;
+rr=0;
+n_star=0;
+y_star=0;
+w_star=0;
+a=0;
+nu=0;
+u=0;
+y_f=0;
+pi_f=0;
+y_e=0;
+pi_e=0;
+pi_exp=0;
+y_exp=0;
+end;
+
+
+steady;
+check;
+
+
+
+ M_.Sigma_e=eye(M_.exo_nbr)*0;
+ shocks;
+ var eps_a; stderr 0.01;
+ end;
+
+
+ varobs yobs piobs robs ;
+
+
+
+stoch_simul(order = 1,irf=40, periods=0, partial_information,contemporaneous_correlation,TeX)y,x,pi,r,rr;
+collect_latex_files;
+
+[status, cmdout]=system(['pdflatex -halt-on-error -interaction=nonstopmode ' M_.fname '_TeX_binder.tex']);
+if status
+ cmdout
+ error('TeX-File did not compile.')
+end