From 36f264a1ca9eaddf455d506d6e6e9a87d2921fdb Mon Sep 17 00:00:00 2001 From: michel Date: Wed, 2 Jul 2008 11:15:19 +0000 Subject: [PATCH] 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 --- matlab/DiffuseLikelihood1_Z.m | 6 +- matlab/DsgeLikelihood.m | 217 +++++++++++++++++-------------- matlab/DsgeSmoother.m | 227 +++++++++++++++++++-------------- matlab/dynare_estimation.m | 6 + matlab/dynare_sensitivity.m | 13 ++ matlab/global_initialization.m | 59 ++++----- matlab/rplot.m | 7 +- matlab/steady.m | 14 +- 8 files changed, 314 insertions(+), 235 deletions(-) diff --git a/matlab/DiffuseLikelihood1_Z.m b/matlab/DiffuseLikelihood1_Z.m index 34cbf2796..525a07428 100644 --- a/matlab/DiffuseLikelihood1_Z.m +++ b/matlab/DiffuseLikelihood1_Z.m @@ -53,8 +53,8 @@ function [LIK, lik] = DiffuseLikelihood1_Z(T,Z,R,Q,Pinf,Pstar,Y,start) return else Fstar = Z*Pstar*Z'; - iFstar = inv(F); - dFstar = det(F); + iFstar = inv(Fstar); + dFstar = det(Fstar); Kstar = Pstar*Z'*iFstar; lik(t) = log(dFstar) + v'*iFstar*v; 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; v = Y(:,t)-Z*a; a = T*(a+K*v); - lik(t) = v*iF*v; + lik(t) = v'*iF*v; end lik(t) = lik(t) + reste*log(dF); diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m index 831f0277b..dae405f89 100644 --- a/matlab/DsgeLikelihood.m +++ b/matlab/DsgeLikelihood.m @@ -151,27 +151,19 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data %------------------------------------------------------------------------------ % 3. Initial condition of the Kalman filter %------------------------------------------------------------------------------ + kalman_algo = options_.kalman_algo; if options_.lik_init == 1 % Kalman filter - Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium); - Pinf = []; + kalman_algo = 1; + Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium); + Pinf = []; elseif options_.lik_init == 2 % Old Diffuse Kalman filter - Pstar = 10*eye(np); - Pinf = []; + kalman_algo = 1; + Pstar = 10*eye(np); + Pinf = []; elseif options_.lik_init == 3 % Diffuse Kalman filter - if options_.kalman_algo < 4 - 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 + if kalman_algo == 1 + kalman_algo == 3 + end [QT,ST] = schur(T); e1 = abs(ordeig(ST)) > 2-options_.qz_criterium; [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); B = QT'*R*Q*R'*QT; for i=np:-1:nk+2 - if ST(i,i-1) == 0 - if i == np - c = zeros(np-nk,1); - else - 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); - end - q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i); - Pstar(nk1:i,i) = q\(B(nk1:i,i)+c); - Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; - else - if i == np - c = zeros(np-nk,1); - c1 = zeros(np-nk,1); - else - 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-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)')+... - 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); - end - 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)]; - 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-1) = z(i-nk+1:end); - Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; - Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)'; - i = i - 1; - end + if ST(i,i-1) == 0 + if i == np + c = zeros(np-nk,1); + else + 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); + end + q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i); + Pstar(nk1:i,i) = q\(B(nk1:i,i)+c); + Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; + else + if i == np + c = zeros(np-nk,1); + c1 = zeros(np-nk,1); + else + 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-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)')+... + 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); + end + 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)]; + 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-1) = z(i-nk+1:end); + Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; + Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)'; + i = i - 1; + end end 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); - Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,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)); end Z = QT(mf,:); 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 %------------------------------------------------------------------------------ % 4. Likelihood evaluation %------------------------------------------------------------------------------ if any(any(H ~= 0)) % should be replaced by a flag - if options_.kalman_algo == 1 - LIK = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,data,trend,start); - if isinf(LIK) & ~estim_params_.ncn %% The univariate approach considered here doesn't - %% apply when H has some off-diagonal elements. - LIK = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,data,trend,start); - elseif isinf(LIK) & estim_params_.ncn - LIK = DiffuseLikelihoodH3corr(T,R,Q,H,Pinf,Pstar,data,trend,start); - end - elseif options_.kalman_algo == 3 - if ~estim_params_.ncn %% The univariate approach considered here doesn't - %% 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 - end + if kalman_algo == 1 + LIK = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,data,trend,start); + if isinf(LIK) + kalman_algo = 2; + if ~estim_params_.ncn + %% The univariate approach considered here doesn't + %% 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 + end + elseif kalman_algo == 2 + if ~estim_params_.ncn + %% The univariate approach considered here doesn't + %% 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 - if options_.kalman_algo == 1 - %nv = size(bayestopt_.Z,1); - %LIK = kalman_filter(bayestopt_.Z,zeros(nv,nv),T,R,Q,data,zeros(size(T,1),1),Pstar,'u'); - %tic - LIK = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,data,trend,start); - %toc - %tic - %LIK1 = Diffuse_Likelihood1(T,R,Q,Pinf,Pstar,data,trend,start); - %toc - %LIK1-LIK - %if abs(LIK1-LIK)>0.0000000001 - % disp(['LIK1 and LIK are not equal! ' num2str(abs(LIK1-LIK))]) - %end - if isinf(LIK) - LIK = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start); - end - elseif options_.kalman_algo == 3 - LIK = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start); - elseif options_.kalman_algo == 4 - 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 + if kalman_algo == 1 + LIK = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,data,trend,start); + if isinf(LIK) + kalman_algo = 2 + LIK = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start); + end + elseif kalman_algo == 2 + LIK = DiffuseLikelihood3(T,R,Q,Pinf,Pstar,data,trend,start); + elseif options_.kalman_algo == 3 + data1 = data - trend; + LIK = DiffuseLikelihood1_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start); + if isinf(LIK) + kalman_algo = 4 + LIK = DiffuseLikelihood3_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start); + end + elseif kalman_algo == 4 + data1 = data - trend; + LIK = DiffuseLikelihood3_Z(ST,Z,R1,Q,Pinf,Pstar,data1,start); + end end if imag(LIK) ~= 0 - likelihood = bayestopt_.penalty; + likelihood = bayestopt_.penalty; else - likelihood = LIK; + likelihood = LIK; end % ------------------------------------------------------------------------------ % Adds prior if necessary % ------------------------------------------------------------------------------ lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4); - fval = (likelihood-lnprior); \ No newline at end of file + fval = (likelihood-lnprior); + options_.kalman_algo = kalman_algo; \ No newline at end of file diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index d5d5b59a7..91c8372ba 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -92,27 +92,17 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d, Q = M_.Sigma_e; H = M_.H; + kalman_algo = options_.kalman_algo; if options_.lik_init == 1 % Kalman filter - Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium); - Pinf = []; + kalman_algo = 1; + Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium); + Pinf = []; elseif options_.lik_init == 2 % Old Diffuse Kalman filter - Pstar = 10*eye(np); - Pinf = []; + kalman_algo = 1; + Pstar = 10*eye(np); + Pinf = []; elseif options_.lik_init == 3 % Diffuse Kalman filter - if options_.kalman_algo < 4 - 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 + kalman_algo = 3; [QT,ST] = schur(T); e1 = abs(ordeig(ST)) > 2-options_.qz_criterium; [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); B = QT'*R*Q*R'*QT; for i=np:-1:nk+2 - if ST(i,i-1) == 0 - if i == np - c = zeros(np-nk,1); - else - 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); - end - q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i); - Pstar(nk1:i,i) = q\(B(nk1:i,i)+c); - Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; - else - if i == np - c = zeros(np-nk,1); - c1 = zeros(np-nk,1); - else - 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-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)')+... - 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); - end - 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)]; - 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-1) = z(i-nk+1:end); - Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; - Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)'; - i = i - 1; - end + if ST(i,i-1) == 0 + if i == np + c = zeros(np-nk,1); + else + 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); + end + q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i); + Pstar(nk1:i,i) = q\(B(nk1:i,i)+c); + Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; + else + if i == np + c = zeros(np-nk,1); + c1 = zeros(np-nk,1); + else + 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-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)')+... + 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); + end + 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)]; + 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-1) = z(i-nk+1:end); + Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; + Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)'; + i = i - 1; + end end 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); - Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,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)); end Z = QT(mf,:); R1 = QT'*R; - end end % ----------------------------------------------------------------------------- % 4. Kalman smoother % ----------------------------------------------------------------------------- - if any(any(H ~= 0)) % should be replaced by a flag - if options_.kalman_algo == 1 - [alphahat,epsilonhat,etahat,ahat,aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf); - if all(alphahat(:)==0) - [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 any(any(H ~= 0)) % should be replaced by a flag + if kalman_algo == 1 + [alphahat,epsilonhat,etahat,ahat,aK] = ... + DiffuseKalmanSmootherH1(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf); if all(alphahat(:)==0) - options_.kalman_algo = 5; - [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ... + 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 + 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, ... data1,nobs,np,smpl); - end - else - [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ... + end + else + [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ... 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 - 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 diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m index ce63b99ec..ed50be843 100644 --- a/matlab/dynare_estimation.m +++ b/matlab/dynare_estimation.m @@ -37,6 +37,12 @@ if options_.order > 1 options_.order = 1; 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 options_.noconstant = 1; end diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m index d49a39fa3..130ec044a 100644 --- a/matlab/dynare_sensitivity.m +++ b/matlab/dynare_sensitivity.m @@ -99,6 +99,19 @@ if ~(exist('stab_map_','file')==6 | exist('stab_map_','file')==2), 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, options_gsa.redform=1; diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 76708d449..a89860f89 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -104,49 +104,50 @@ function global_initialization() options_.ramsey_policy = 0; % 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_.prefilter = 0; - options_.presample = 0; + options_.kalman_algo = 0; + options_.kalman_tol = 1e-12; options_.lik_algo = 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_jscale = 0.2; 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_nblck = 2; - options_.load_mh_file = 0; options_.mh_recover = 0; - options_.nodiagnostic = 0; - options_.loglinear = 0; - options_.unit_root_vars = []; - options_.bayesian_irf = 0; - options_.bayesian_th_moments = 0; - options_.smoother = 0; + options_.mh_replic = 20000; + options_.mode_check = 0; + options_.mode_compute = 4; + options_.mode_file = ''; 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_.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_.cutoff = 1e-12; + options_.smoother = 0; options_.student_degrees_of_freedom = 3; options_.subdraws = []; - options_.PosteriorSampleSize = 1000; - options_.MaximumNumberOfMegaBytes = 111; + options_.unit_root_vars = []; + options_.use_mh_covariance_matrix = 0; % Misc options_.conf_sig = 0.6; diff --git a/matlab/rplot.m b/matlab/rplot.m index b35b636af..b0f0d2959 100644 --- a/matlab/rplot.m +++ b/matlab/rplot.m @@ -35,10 +35,10 @@ for k=1:size(s1,1) y = [y; oo_.endo_simul(strmatch(s1(k,:),M_.endo_names,'exact'),:)] ; end -if options_.smpl == 0 +if options_.dsample == 0 i = [M_.maximum_lag:size(oo_.endo_simul,2)]' ; 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 t = ['Plot of '] ; @@ -75,9 +75,6 @@ elseif rplottype == 2 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 diff --git a/matlab/steady.m b/matlab/steady.m index 7d3950af1..bbe078c00 100644 --- a/matlab/steady.m +++ b/matlab/steady.m @@ -37,7 +37,19 @@ function steady() homotopy3(options_.homotopy_values, options_.homotopy_steps); 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('STEADY-STATE RESULTS:')