stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

put back DsgeLikelihood_hh.m that is used by option 5 of the optimizer 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1918 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
michel 2008-07-02 11:15:19 +00:00
parent 9074fdba91
commit 36f264a1ca
8 changed files with 314 additions and 235 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
matlab/DiffuseLikelihood1_Z.m
View File

@ -53,8 +53,8 @@ function [LIK, lik] = DiffuseLikelihood1_Z(T,Z,R,Q,Pinf,Pstar,Y,start)
return return
else else
Fstar = Z*Pstar*Z'; Fstar = Z*Pstar*Z';
iFstar = inv(F); iFstar = inv(Fstar);
dFstar = det(F); dFstar = det(Fstar);
Kstar = Pstar*Z'*iFstar; Kstar = Pstar*Z'*iFstar;
lik(t) = log(dFstar) + v'*iFstar*v; lik(t) = log(dFstar) + v'*iFstar*v;
Pinf = T*Pinf*transpose(T); Pinf = T*Pinf*transpose(T);
@ -109,7 +109,7 @@ function [LIK, lik] = DiffuseLikelihood1_Z(T,Z,R,Q,Pinf,Pstar,Y,start)
t = t+1; t = t+1;
v = Y(:,t)-Z*a; v = Y(:,t)-Z*a;
a = T*(a+K*v); a = T*(a+K*v);
lik(t) = v*iF*v; lik(t) = v'*iF*v;
end end
lik(t) = lik(t) + reste*log(dF); lik(t) = lik(t) + reste*log(dF);

217
matlab/DsgeLikelihood.m
View File

@ -151,27 +151,19 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter % 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
kalman_algo = options_.kalman_algo;
if options_.lik_init == 1 % Kalman filter if options_.lik_init == 1 % Kalman filter
Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium); kalman_algo = 1;
Pinf = []; Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium);
Pinf = [];
elseif options_.lik_init == 2 % Old Diffuse Kalman filter elseif options_.lik_init == 2 % Old Diffuse Kalman filter
Pstar = 10*eye(np); kalman_algo = 1;
Pinf = []; Pstar = 10*eye(np);
Pinf = [];
elseif options_.lik_init == 3 % Diffuse Kalman filter elseif options_.lik_init == 3 % Diffuse Kalman filter
if options_.kalman_algo < 4 if kalman_algo == 1
Pstar = zeros(np,np); kalman_algo == 3
ivs = bayestopt_.restrict_var_list_stationary; end
R1 = R(ivs,:);
Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R1*Q*R1',options_.qz_criterium);
% Pinf = bayestopt_.Pinf;
% by M. Ratto
RR=T(:,bayestopt_.restrict_var_list_nonstationary);
i=find(abs(RR)>1.e-10);
R0=zeros(size(RR));
R0(i)=sign(RR(i));
Pinf=R0*R0';
% by M. Ratto
else
[QT,ST] = schur(T); [QT,ST] = schur(T);
e1 = abs(ordeig(ST)) > 2-options_.qz_criterium; e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
[QT,ST] = ordschur(QT,ST,e1); [QT,ST] = ordschur(QT,ST,e1);
@ -183,104 +175,131 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
Pstar = zeros(np,np); Pstar = zeros(np,np);
B = QT'*R*Q*R'*QT; B = QT'*R*Q*R'*QT;
for i=np:-1:nk+2 for i=np:-1:nk+2
if ST(i,i-1) == 0 if ST(i,i-1) == 0
if i == np if i == np
c = zeros(np-nk,1); c = zeros(np-nk,1);
else else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+... c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i); ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end end
q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i); q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
Pstar(nk1:i,i) = q\(B(nk1:i,i)+c); Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
else else
if i == np if i == np
c = zeros(np-nk,1); c = zeros(np-nk,1);
c1 = zeros(np-nk,1); c1 = zeros(np-nk,1);
else else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+... c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+... ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1); ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+... c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+... ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i); ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end end
q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);... q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
-ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)]; -ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1]; z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
Pstar(nk1:i,i) = z(1:(i-nk)); Pstar(nk1:i,i) = z(1:(i-nk));
Pstar(nk1:i,i-1) = z(i-nk+1:end); Pstar(nk1:i,i-1) = z(i-nk+1:end);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)'; Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
i = i - 1; i = i - 1;
end end
end end
if i == nk+2 if i == nk+2
c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1); c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1)); Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
end end
Z = QT(mf,:); Z = QT(mf,:);
R1 = QT'*R; R1 = QT'*R;
end [u,s,v]=svd(Z*ST(:,1:nk),0);
k = find(abs(diag(s)) < 1e-8);
if length(k) > 0
[junk,k1] = max(abs(v(:,k)));
dd =ones(nk,1);
dd(k1) = zeros(length(k1),1);
Pinf(1:nk,1:nk) = diag(dd);
end
end end
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
% 4. Likelihood evaluation % 4. Likelihood evaluation
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
if any(any(H ~= 0)) % should be replaced by a flag if any(any(H ~= 0)) % should be replaced by a flag
if options_.kalman_algo == 1 if kalman_algo == 1
LIK = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,data,trend,start); LIK = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,data,trend,start);
if isinf(LIK) & ~estim_params_.ncn %% The univariate approach considered here doesn't if isinf(LIK)
%% apply when H has some off-diagonal elements. kalman_algo = 2;
LIK = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,data,trend,start); if ~estim_params_.ncn
elseif isinf(LIK) & estim_params_.ncn %% The univariate approach considered here doesn't
LIK = DiffuseLikelihoodH3corr(T,R,Q,H,Pinf,Pstar,data,trend,start); %% apply when H has some off-diagonal elements.
end LIK = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,data,trend,start);
elseif options_.kalman_algo == 3 else
if ~estim_params_.ncn %% The univariate approach considered here doesn't LIK = DiffuseLikelihoodH3corr(T,R,Q,H,Pinf,Pstar,data,trend,start);
%% apply when H has some off-diagonal elements. end
LIK = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,data,trend,start); end
else elseif kalman_algo == 2
LIK = DiffuseLikelihoodH3corr(T,R,Q,H,Pinf,Pstar,data,trend,start); if ~estim_params_.ncn
end %% The univariate approach considered here doesn't
end %% apply when H has some off-diagonal elements.
LIK = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,data,trend,start);
else
LIK = DiffuseLikelihoodH3corr(T,R,Q,H,Pinf,Pstar,data,trend,start);
end
elseif kalman_algo == 3
data1 = data - trend;
LIK = DiffuseLikelihoodH1_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,trend,start);
if isinf(LIK)
kalman_algo = 4;
if ~estim_params_.ncn
%% The univariate approach considered here doesn't
%% apply when H has some off-diagonal elements.
LIK = DiffuseLikelihoodH3_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,trend,start);
else
LIK = DiffuseLikelihoodH3corr_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,trend,start);
end
end
elseif kalman_algo == 4
data1 = data - trend;
if ~estim_params_.ncn
%% The univariate approach considered here doesn't
%% apply when H has some off-diagonal elements.
LIK = DiffuseLikelihoodH3_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,trend,start);
else
LIK = DiffuseLikelihoodH3corr_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,trend,start);
end
end
else else
if options_.kalman_algo == 1 if kalman_algo == 1
%nv = size(bayestopt_.Z,1); LIK = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,data,trend,start);
%LIK = kalman_filter(bayestopt_.Z,zeros(nv,nv),T,R,Q,data,zeros(size(T,1),1),Pstar,'u'); if isinf(LIK)
%tic kalman_algo = 2
LIK = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,data,trend,start); LIK = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start);
%toc end
%tic elseif kalman_algo == 2
%LIK1 = Diffuse_Likelihood1(T,R,Q,Pinf,Pstar,data,trend,start); LIK = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start);
%toc elseif options_.kalman_algo == 3
%LIK1-LIK data1 = data - trend;
%if abs(LIK1-LIK)>0.0000000001 LIK = DiffuseLikelihood1_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start);
% disp(['LIK1 and LIK are not equal! ' num2str(abs(LIK1-LIK))]) if isinf(LIK)
%end kalman_algo = 4
if isinf(LIK) LIK = DiffuseLikelihood3_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start);
LIK = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start); end
end elseif kalman_algo == 4
elseif options_.kalman_algo == 3 data1 = data - trend;
LIK = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start); LIK = DiffuseLikelihood3_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start);
elseif options_.kalman_algo == 4 end
data1 = data - trend;
LIK = DiffuseLikelihood1_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start);
if isinf(LIK)
LIK = DiffuseLikelihood3_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start);
end
elseif options_.kalman_algo == 5
data1 = data - trend;
LIK = DiffuseLikelihood3_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start);
end
end end
if imag(LIK) ~= 0 if imag(LIK) ~= 0
likelihood = bayestopt_.penalty; likelihood = bayestopt_.penalty;
else else
likelihood = LIK; likelihood = LIK;
end end
% ------------------------------------------------------------------------------ % ------------------------------------------------------------------------------
% Adds prior if necessary % Adds prior if necessary
% ------------------------------------------------------------------------------ % ------------------------------------------------------------------------------
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4); lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
fval = (likelihood-lnprior); fval = (likelihood-lnprior);
options_.kalman_algo = kalman_algo;

227
matlab/DsgeSmoother.m
View File

@ -92,27 +92,17 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
Q = M_.Sigma_e; Q = M_.Sigma_e;
H = M_.H; H = M_.H;
kalman_algo = options_.kalman_algo;
if options_.lik_init == 1 % Kalman filter if options_.lik_init == 1 % Kalman filter
Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium); kalman_algo = 1;
Pinf = []; Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium);
Pinf = [];
elseif options_.lik_init == 2 % Old Diffuse Kalman filter elseif options_.lik_init == 2 % Old Diffuse Kalman filter
Pstar = 10*eye(np); kalman_algo = 1;
Pinf = []; Pstar = 10*eye(np);
Pinf = [];
elseif options_.lik_init == 3 % Diffuse Kalman filter elseif options_.lik_init == 3 % Diffuse Kalman filter
if options_.kalman_algo < 4 kalman_algo = 3;
Pstar = zeros(np,np);
ivs = bayestopt_.restrict_var_list_stationary;
R1 = R(ivs,:);
Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R1*Q*R1',options_.qz_criterium);
% Pinf = bayestopt_.Pinf;
% by M. Ratto
RR=T(:,bayestopt_.restrict_var_list_nonstationary);
i=find(abs(RR)>1.e-10);
R0=zeros(size(RR));
R0(i)=sign(RR(i));
Pinf=R0*R0';
% by M. Ratto
else
[QT,ST] = schur(T); [QT,ST] = schur(T);
e1 = abs(ordeig(ST)) > 2-options_.qz_criterium; e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
[QT,ST] = ordschur(QT,ST,e1); [QT,ST] = ordschur(QT,ST,e1);
@ -124,99 +114,140 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
Pstar = zeros(np,np); Pstar = zeros(np,np);
B = QT'*R*Q*R'*QT; B = QT'*R*Q*R'*QT;
for i=np:-1:nk+2 for i=np:-1:nk+2
if ST(i,i-1) == 0 if ST(i,i-1) == 0
if i == np if i == np
c = zeros(np-nk,1); c = zeros(np-nk,1);
else else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+... c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i); ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end end
q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i); q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
Pstar(nk1:i,i) = q\(B(nk1:i,i)+c); Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
else else
if i == np if i == np
c = zeros(np-nk,1); c = zeros(np-nk,1);
c1 = zeros(np-nk,1); c1 = zeros(np-nk,1);
else else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+... c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+... ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1); ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+... c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+... ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i); ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end end
q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);... q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
-ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)]; -ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1]; z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
Pstar(nk1:i,i) = z(1:(i-nk)); Pstar(nk1:i,i) = z(1:(i-nk));
Pstar(nk1:i,i-1) = z(i-nk+1:end); Pstar(nk1:i,i-1) = z(i-nk+1:end);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)'; Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
i = i - 1; i = i - 1;
end end
end end
if i == nk+2 if i == nk+2
c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1); c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1)); Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
end end
Z = QT(mf,:); Z = QT(mf,:);
R1 = QT'*R; R1 = QT'*R;
end
end end
% ----------------------------------------------------------------------------- % -----------------------------------------------------------------------------
% 4. Kalman smoother % 4. Kalman smoother
% ----------------------------------------------------------------------------- % -----------------------------------------------------------------------------
if any(any(H ~= 0)) % should be replaced by a flag if any(any(H ~= 0)) % should be replaced by a flag
if options_.kalman_algo == 1 if kalman_algo == 1
[alphahat,epsilonhat,etahat,ahat,aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf); [alphahat,epsilonhat,etahat,ahat,aK] = ...
if all(alphahat(:)==0) DiffuseKalmanSmootherH1(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
[alphahat,epsilonhat,etahat,ahat,aK] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
end
elseif options_.kalman_algo == 3
[alphahat,epsilonhat,etahat,ahat,aK] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
end
else
if options_.kalman_algo == 1
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
if all(alphahat(:)==0)
options_.kalman_algo = 3;
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
end
elseif options_.kalman_algo == 3
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
elseif options_.kalman_algo == 4 | options_.kalman_algo == 5
data1 = Y - trend;
if options_.kalman_algo == 4
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother1_Z(ST, ...
Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl);
if all(alphahat(:)==0) if all(alphahat(:)==0)
options_.kalman_algo = 5; kalman_algo = 2;
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ... if ~estim_params.ncn
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
else
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf,Pstar,Y,trend, ...
nobs,np,smpl,mf);
end
end
elseif options_.kalman_algo == 2
if ~estim_params.ncn
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
else
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf,Pstar,Y,trend, ...
nobs,np,smpl,mf);
end
elseif kalman_algo == 3
data1 = Y - trend;
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH1_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,nobs,np,smpl);
if all(alphahat(:)==0)
kalman_algo = 4;
if ~estim_params.ncn
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH3_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1,nobs,np,smpl);
else
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH3corr_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1, ...
nobs,np,smpl);
end
end
elseif options_.kalman_algo == 4
data1 = Y - trend;
if ~estim_params.ncn
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH3_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1, ...
nobs,np,smpl);
else
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH3corr_Z(ST,Z,R1,Q,H,Pinf,Pstar,data1, ...
nobs,np,smpl);
end
end
else
if kalman_algo == 1
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
if all(alphahat(:)==0)
kalman_algo = 2;
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
end
elseif kalman_algo == 2
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
elseif kalman_algo == 3 | kalman_algo == 4
data1 = Y - trend;
if options_.kalman_algo == 3
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother1_Z(ST, ...
Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl);
if all(alphahat(:)==0)
options_.kalman_algo = 4;
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ...
Z,R1,Q,Pinf,Pstar, ... Z,R1,Q,Pinf,Pstar, ...
data1,nobs,np,smpl); data1,nobs,np,smpl);
end end
else else
[alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ... [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ...
Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl); Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl);
end
alphahat = QT*alphahat;
ahat = QT*ahat;
if options_.nk > 0
nk = options_.nk;
for jnk=1:nk
aK(jnk,:,:) = QT*squeeze(aK(jnk,:,:));
for i=1:size(PK,4)
PK(jnk,:,:,i) = QT*squeeze(PK(jnk,:,:,i))*QT';
end
for i=1:size(decomp,4)
decomp(jnk,:,:,i) = QT*squeeze(decomp(jnk,:,:,i));
end
end
for i=1:size(P,4)
P(:,:,i) = QT*squeeze(P(:,:,i))*QT';
end
end
end end
alphahat = QT*alphahat;
ahat = QT*ahat;
if options_.nk > 0
nk = options_.nk;
for jnk=1:nk
aK(jnk,:,:) = QT*squeeze(aK(jnk,:,:));
for i=1:size(PK,4)
PK(jnk,:,:,i) = QT*squeeze(PK(jnk,:,:,i))*QT';
end
for i=1:size(decomp,4)
decomp(jnk,:,:,i) = QT*squeeze(decomp(jnk,:,:,i));
end
end
for i=1:size(P,4)
P(:,:,i) = QT*squeeze(P(:,:,i))*QT';
end
end
end
end end

6
matlab/dynare_estimation.m
View File

@ -37,6 +37,12 @@ if options_.order > 1
options_.order = 1; options_.order = 1;
end end
if (~isempty(options_.unit_root_vars) || options_.diffuse_filter == 1)
if options_.lik_init == 1
options_.lik_init = 3
end
end
if options_.prefilter == 1 if options_.prefilter == 1
options_.noconstant = 1; options_.noconstant = 1;
end end

13
matlab/dynare_sensitivity.m
View File

@ -99,6 +99,19 @@ if ~(exist('stab_map_','file')==6 | exist('stab_map_','file')==2),
end end
end end
if ~(exist('gsa_sdp','file')==6 | exist('gsa_sdp','file')==2),
dynare_root = strrep(which('dynare.m'),'dynare.m','');
ss_anova_path = [dynare_root 'gsa/ss_anova_recurs'];
if exist(ss_anova_path)
addpath(ss_anova_path,path)
else
disp('Download Dynare sensitivity routines at:')
disp('http://eemc.jrc.ec.europa.eu/softwareDYNARE-Dowload.htm')
disp(' ' )
error('Mapping routines missing!')
end
end
if options_gsa.morris, if options_gsa.morris,
options_gsa.redform=1; options_gsa.redform=1;

59
matlab/global_initialization.m
View File

@ -104,49 +104,50 @@ function global_initialization()
options_.ramsey_policy = 0; options_.ramsey_policy = 0;
% estimation % estimation
options_.load_mh_file = 0; options_.MaxNumberOfBytes = 1e6;
options_.MaximumNumberOfMegaBytes = 111;
options_.PosteriorSampleSize = 1000;
options_.bayesian_irf = 0;
options_.bayesian_th_moments = 0;
options_.cutoff = 1e-12;
options_.diffuse_d = [];
options_.diffuse_filter = 0;
options_.filter_step_ahead = [];
options_.filtered_vars = 0;
options_.first_obs = 1; options_.first_obs = 1;
options_.prefilter = 0; options_.kalman_algo = 0;
options_.presample = 0; options_.kalman_tol = 1e-12;
options_.lik_algo = 1; options_.lik_algo = 1;
options_.lik_init = 1; options_.lik_init = 1;
options_.mh_replic = 20000; options_.load_mh_file = 0;
options_.load_mh_file = 0;
options_.logdata = 0;
options_.loglinear = 0;
options_.markowitz = 0.5;
options_.mh_conf_sig = 0.90;
options_.mh_drop = 0.5; options_.mh_drop = 0.5;
options_.mh_jscale = 0.2; options_.mh_jscale = 0.2;
options_.mh_init_scale = 2*options_.mh_jscale; options_.mh_init_scale = 2*options_.mh_jscale;
options_.mode_file = '';
options_.mode_compute = 4;
options_.mode_check = 0;
options_.prior_trunc = 1e-10;
options_.mh_conf_sig = 0.90;
options_.mh_mode = 1; options_.mh_mode = 1;
options_.mh_nblck = 2; options_.mh_nblck = 2;
options_.load_mh_file = 0;
options_.mh_recover = 0; options_.mh_recover = 0;
options_.nodiagnostic = 0; options_.mh_replic = 20000;
options_.loglinear = 0; options_.mode_check = 0;
options_.unit_root_vars = []; options_.mode_compute = 4;
options_.bayesian_irf = 0; options_.mode_file = '';
options_.bayesian_th_moments = 0;
options_.smoother = 0;
options_.moments_varendo = 0; options_.moments_varendo = 0;
options_.filtered_vars = 0;
options_.kalman_algo = 1;
options_.kalman_tol = 1e-12;
options_.posterior_mode_estimation = 1;
options_.MaxNumberOfBytes = 1e6;
options_.filter_step_ahead = [];
options_.diffuse_d = [];
options_.logdata = 0;
options_.use_mh_covariance_matrix = 0;
options_.noconstant = 0; options_.noconstant = 0;
options_.markowitz = 0.5; options_.nodiagnostic = 0;
options_.posterior_mode_estimation = 1;
options_.prefilter = 0;
options_.presample = 0;
options_.prior_trunc = 1e-10;
options_.simulation_method = 0; options_.simulation_method = 0;
options_.cutoff = 1e-12; options_.smoother = 0;
options_.student_degrees_of_freedom = 3; options_.student_degrees_of_freedom = 3;
options_.subdraws = []; options_.subdraws = [];
options_.PosteriorSampleSize = 1000; options_.unit_root_vars = [];
options_.MaximumNumberOfMegaBytes = 111; options_.use_mh_covariance_matrix = 0;
% Misc % Misc
options_.conf_sig = 0.6; options_.conf_sig = 0.6;

7
matlab/rplot.m
View File

@ -35,10 +35,10 @@ for k=1:size(s1,1)
y = [y; oo_.endo_simul(strmatch(s1(k,:),M_.endo_names,'exact'),:)] ; y = [y; oo_.endo_simul(strmatch(s1(k,:),M_.endo_names,'exact'),:)] ;
end end
if options_.smpl == 0 if options_.dsample == 0
i = [M_.maximum_lag:size(oo_.endo_simul,2)]' ; i = [M_.maximum_lag:size(oo_.endo_simul,2)]' ;
else else
i = [options_.smpl(1)+M_.maximum_lag:options_.smpl(2)+M_.maximum_lag]' ; i = [options_.dsample(1)+M_.maximum_lag:options_.dsample(2)+M_.maximum_lag]' ;
end end
t = ['Plot of '] ; t = ['Plot of '] ;
@ -75,9 +75,6 @@ elseif rplottype == 2
end end
end end
% 02/28/01 MJ replaced bseastr by MATLAB's strmatch
% 06/19/01 MJ added 'exact' to strmatch calls
% 06/25/03 MJ correction when options_.smpl ~= 0

14
matlab/steady.m
View File

@ -37,7 +37,19 @@ function steady()
homotopy3(options_.homotopy_values, options_.homotopy_steps); homotopy3(options_.homotopy_values, options_.homotopy_steps);
end end
steady_; if options_.homotopy_mode == 1
homotopy1(options_.homotopy_param,options_.homotopy_exo, ...
options_.homotopy_exodet,options_.homotopy_steps);
elseif options_.homotopy_mode == 2
homotopy2(options_.homotopy_param,options_.homotopy_exo, ...
options_.homotopy_exodet,options_.homotopy_steps);
elseif options_.homotopy_mode == 3
homotopy3(options_.homotopy_param,options_.homotopy_exo, ...
options_.homotopy_exodet,options_.homotopy_steps);
else
steady_;
end
disp(' ') disp(' ')
disp('STEADY-STATE RESULTS:') disp('STEADY-STATE RESULTS:')