change file format to unix
Houtan Bastani
parent
f44db9370d
commit
e4546ba32f
|
|
@ -1,197 +1,197 @@
|
|||
function ms_mardd(options_)
|
||||
% Applies to both linear and exclusion restrictions.
|
||||
% (1) Marginal likelihood function p(Y) for constant structural VAR models, using Chib (1995)'s ``Marginal Likelihood from the Gibbs Output'' in JASA.
|
||||
% (2) Conditional likelihood function f(Y|A0, A+) on the ML estimate for constant exclusion-identified models.
|
||||
% See Forecast (II) pp.67-80.
|
||||
%
|
||||
% Tao Zha, September 1999. Quick revisions, May 2003. Final revision, September 2004.
|
||||
|
||||
% Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
msstart2 % start the program in which everyhting is initialized through msstart2.m
|
||||
if ~options_.ms.indxestima
|
||||
warning(' ')
|
||||
disp('You must set IxEstima=1 in msstart to run this program')
|
||||
disp('Press ctrl-c to abort now')
|
||||
pause
|
||||
end
|
||||
|
||||
A0xhat = zeros(size(A0hat));
|
||||
Apxhat = zeros(size(Aphat));
|
||||
if (0)
|
||||
%Robustness check to see if the same result is obtained with the purterbation of the parameters.
|
||||
for k=1:nvar
|
||||
bk = Uiconst{k}'*A0hat(:,k);
|
||||
gk = Viconst{k}'*Aphat(:,k);
|
||||
A0xhat(:,k) = Uiconst{k}*(bk + 5.2*randn(size(bk))); % Perturbing the posterior estimate.
|
||||
Apxhat(:,k) = Viconst{k}*(gk + 5.2*randn(size(gk))); % Perturbing the posterior estimate.
|
||||
end
|
||||
else
|
||||
%At the posterior estimate.
|
||||
A0xhat = A0hat; % ML estimate of A0
|
||||
Apxhat = Aphat; % ML estimate of A+
|
||||
end
|
||||
%--- Rename variables.
|
||||
YatYa = yty;
|
||||
XatYa = xty;
|
||||
ytx = xty';
|
||||
YatXa = ytx;
|
||||
XatXa = xtx;
|
||||
|
||||
|
||||
|
||||
%--------- The log value of p(A0,A+) at some point such as the peak ----------
|
||||
vlog_a0p = 0;
|
||||
Yexpt=0; % exponential term for Y in p(Y|A0,A+) at some point such as the peak
|
||||
Apexpt=0.0; % 0.0 because we have chosen posterior estimate of A+ as A+*. Exponential term for A+ conditional on A0 and Y
|
||||
%======= Computing the log prior pdf of a0a+ and the exponential term for Y in p(Y|A0,A+).
|
||||
for k=1:nvar
|
||||
a0k = A0xhat(:,k); % meaningful parameters in the kth equation.
|
||||
apk = Apxhat(:,k); % meaningful parameters in the kth equation.
|
||||
|
||||
%--- Prior settings.
|
||||
S0bar = H0invtld{k}; %See Claim 2 on p.69b.
|
||||
Spbar = Hpinvtld{k};
|
||||
bk = Uiconst{k}'*a0k; % free parameters in the kth equation.
|
||||
gk = Viconst{k}'*apk; % free parameters in the kth equation.
|
||||
gbark = Ptld{k}*bk; % bar: prior
|
||||
|
||||
%--- The exponential term for Y in p(Y|A0,A+)
|
||||
Yexpt = Yexpt - 0.5*(a0k'*YatYa*a0k - 2*apk'*XatYa*a0k + apk'*XatXa*apk);
|
||||
%--- The log prior pdf.
|
||||
vlog_a0p = vlog_a0p - 0.5*(size(Uiconst{k},2)+size(Viconst{k},2))*log(2*pi) + 0.5*log(abs(det(S0bar))) + ...
|
||||
0.5*log(abs(det(Spbar))) - 0.5*(bk'*S0bar*bk+(gk-gbark)'*Spbar*(gk-gbark));
|
||||
%--- For p(A+|Y,a0) only.
|
||||
tmpd = gk - Pmat{k}*bk;
|
||||
Apexpt = Apexpt - 0.5*tmpd'*(Hpinv{k}*tmpd);
|
||||
end
|
||||
vlog_a0p
|
||||
|
||||
%--------- The log value of p(Y|A0,A+) at some point such as the peak. ----------
|
||||
%--------- Note that logMarLHres is the same as vlog_Y_a, just to double check. ----------
|
||||
vlog_Y_a = -0.5*nvar*fss*log(2*pi) + fss*log(abs(det(A0xhat))) + Yexpt
|
||||
% a: given a0 and a+
|
||||
logMarLHres = 0; % Initialize log of the marginal likelihood (restricted or constant parameters).
|
||||
for ki=1:fss %ndobs+1:fss % Forward recursion to get the marginal likelihood. See F on p.19 and pp. 48-49.
|
||||
%---- Restricted log marginal likelihood function (constant parameters).
|
||||
[A0l,A0u] = lu(A0xhat);
|
||||
ada = sum(log(abs(diag(A0u)))); % log|A0|
|
||||
termexp = y(ki,:)*A0xhat - phi(ki,:)*Apxhat; % 1-by-nvar
|
||||
logMarLHres = logMarLHres - (0.5*nvar)*log(2*pi) + ada - 0.5*termexp*termexp'; % log value
|
||||
end
|
||||
logMarLHres
|
||||
|
||||
|
||||
%--------- The log value of p(A+|Y,A0) at some point such as the peak ----------
|
||||
totparsp = 0.0;
|
||||
tmpd = 0.0;
|
||||
for k=1:nvar
|
||||
totparsp = totparsp + size(Viconst{k},2);
|
||||
tmpd = tmpd + 0.5*log(abs(det(Hpinv{k})));
|
||||
end
|
||||
vlog_ap_Ya0 = -0.5*totparsp*log(2*pi) + tmpd + Apexpt;
|
||||
|
||||
|
||||
|
||||
|
||||
%===================================
|
||||
% Compute p(a0,k|Y,ao) at some point such as the peak (in this situation, we simply
|
||||
% generate results from the original Gibbs sampler). See FORECAST (2) pp.70-71
|
||||
%===================================
|
||||
%--- Global set up for Gibbs.
|
||||
[Tinv,UT] = fn_gibbsrvar_setup(H0inv, Uiconst, Hpinv, Pmat, Viconst, nvar, fss);
|
||||
%
|
||||
vlog_a0_Yao = zeros(nvar,1);
|
||||
% the log value of p(a0k|Y,ao) where ao: other a's at some point such as the peak of ONLY some a0's
|
||||
vlog=zeros(ndraws2,1);
|
||||
for k=1:nvar
|
||||
bk = Uiconst{k}'*A0xhat(:,k);
|
||||
indx_ks=[k:nvar]; % the columns that exclude 1-(k-1)th columns
|
||||
A0gbs0 = A0hat; % starting at some point such as the peak
|
||||
nk = n0(k);
|
||||
|
||||
if k<nvar
|
||||
%--------- The 1st set of draws to be tossed away. ------------------
|
||||
for draws = 1:ndraws1
|
||||
if ~mod(draws,nbuffer)
|
||||
disp(' ')
|
||||
disp(sprintf('The %dth column or equation in A0 with %d 1st tossed-away draws in Gibbs',k,draws))
|
||||
end
|
||||
A0gbs1 = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
|
||||
A0gbs0=A0gbs1; % repeat the Gibbs sampling
|
||||
end
|
||||
|
||||
|
||||
%--------- The 2nd set of draws to be used. ------------------
|
||||
for draws = 1:ndraws2
|
||||
if ~mod(draws,nbuffer)
|
||||
disp(' ')
|
||||
disp(sprintf('The %dth column or equation in A0 with %d usable draws in Gibbs',k,draws))
|
||||
end
|
||||
[A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
|
||||
%------ See p.71, Forecast (II).
|
||||
%------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
|
||||
Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II).
|
||||
gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
|
||||
[Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
|
||||
%
|
||||
vlog(draws) = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
|
||||
0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
|
||||
0.5*fss*bk'*(Vtr\(Vtr'\bk));
|
||||
|
||||
A0gbs0=A0gbs1; % repeat the Gibbs sampling
|
||||
end
|
||||
vlogm=max(vlog);
|
||||
qlog=vlog-vlogm;
|
||||
vlogxhat=vlogm-log(ndraws2)+log(sum(exp(qlog)));
|
||||
vlog_a0_Yao(k) = vlogxhat;
|
||||
% The log value of p(a0_k|Y,a_others) where a_others: other a's at some point such as the peak of ONLY some a0's
|
||||
else
|
||||
disp(' ')
|
||||
disp(sprintf('The last(6th) column or equation in A0 with no Gibbs draws'))
|
||||
[A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks)
|
||||
%------ See p.71, Forecast (II).
|
||||
%------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
|
||||
Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II).
|
||||
gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
|
||||
[Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
|
||||
%
|
||||
vloglast = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
|
||||
0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
|
||||
0.5*fss*bk'*(Vtr\(Vtr'\bk));
|
||||
vlog_a0_Yao(k) = vloglast;
|
||||
end
|
||||
end
|
||||
ndraws2
|
||||
|
||||
disp('Prior pdf -- log(p(a0hat, a+hat)):');
|
||||
vlog_a0p
|
||||
disp('LH pdf -- log(p(Y|a0hat, a+hat)):');
|
||||
vlog_Y_a
|
||||
disp('Posterior Kernal -- logp(ahat) + logp(Y|ahat):');
|
||||
vlog_Y_a + vlog_a0p
|
||||
disp('Posterior pdf -- log(p(a0_i_hat|a0_other_hat, Y)):');
|
||||
vlog_a0_Yao
|
||||
disp('Posterior pdf -- log(p(aphat|a0hat, Y)):');
|
||||
vlog_ap_Ya0
|
||||
|
||||
%--------- The value of marginal density p(Y) ----------
|
||||
disp(' ');
|
||||
disp(' ');
|
||||
disp('************ Marginal Likelihood of Y or Marginal Data Density: ************');
|
||||
vlogY = vlog_a0p+vlog_Y_a-sum(vlog_a0_Yao)-vlog_ap_Ya0
|
||||
function ms_mardd(options_)
|
||||
% Applies to both linear and exclusion restrictions.
|
||||
% (1) Marginal likelihood function p(Y) for constant structural VAR models, using Chib (1995)'s ``Marginal Likelihood from the Gibbs Output'' in JASA.
|
||||
% (2) Conditional likelihood function f(Y|A0, A+) on the ML estimate for constant exclusion-identified models.
|
||||
% See Forecast (II) pp.67-80.
|
||||
%
|
||||
% Tao Zha, September 1999. Quick revisions, May 2003. Final revision, September 2004.
|
||||
|
||||
% Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
msstart2 % start the program in which everyhting is initialized through msstart2.m
|
||||
if ~options_.ms.indxestima
|
||||
warning(' ')
|
||||
disp('You must set IxEstima=1 in msstart to run this program')
|
||||
disp('Press ctrl-c to abort now')
|
||||
pause
|
||||
end
|
||||
|
||||
A0xhat = zeros(size(A0hat));
|
||||
Apxhat = zeros(size(Aphat));
|
||||
if (0)
|
||||
%Robustness check to see if the same result is obtained with the purterbation of the parameters.
|
||||
for k=1:nvar
|
||||
bk = Uiconst{k}'*A0hat(:,k);
|
||||
gk = Viconst{k}'*Aphat(:,k);
|
||||
A0xhat(:,k) = Uiconst{k}*(bk + 5.2*randn(size(bk))); % Perturbing the posterior estimate.
|
||||
Apxhat(:,k) = Viconst{k}*(gk + 5.2*randn(size(gk))); % Perturbing the posterior estimate.
|
||||
end
|
||||
else
|
||||
%At the posterior estimate.
|
||||
A0xhat = A0hat; % ML estimate of A0
|
||||
Apxhat = Aphat; % ML estimate of A+
|
||||
end
|
||||
%--- Rename variables.
|
||||
YatYa = yty;
|
||||
XatYa = xty;
|
||||
ytx = xty';
|
||||
YatXa = ytx;
|
||||
XatXa = xtx;
|
||||
|
||||
|
||||
|
||||
%--------- The log value of p(A0,A+) at some point such as the peak ----------
|
||||
vlog_a0p = 0;
|
||||
Yexpt=0; % exponential term for Y in p(Y|A0,A+) at some point such as the peak
|
||||
Apexpt=0.0; % 0.0 because we have chosen posterior estimate of A+ as A+*. Exponential term for A+ conditional on A0 and Y
|
||||
%======= Computing the log prior pdf of a0a+ and the exponential term for Y in p(Y|A0,A+).
|
||||
for k=1:nvar
|
||||
a0k = A0xhat(:,k); % meaningful parameters in the kth equation.
|
||||
apk = Apxhat(:,k); % meaningful parameters in the kth equation.
|
||||
|
||||
%--- Prior settings.
|
||||
S0bar = H0invtld{k}; %See Claim 2 on p.69b.
|
||||
Spbar = Hpinvtld{k};
|
||||
bk = Uiconst{k}'*a0k; % free parameters in the kth equation.
|
||||
gk = Viconst{k}'*apk; % free parameters in the kth equation.
|
||||
gbark = Ptld{k}*bk; % bar: prior
|
||||
|
||||
%--- The exponential term for Y in p(Y|A0,A+)
|
||||
Yexpt = Yexpt - 0.5*(a0k'*YatYa*a0k - 2*apk'*XatYa*a0k + apk'*XatXa*apk);
|
||||
%--- The log prior pdf.
|
||||
vlog_a0p = vlog_a0p - 0.5*(size(Uiconst{k},2)+size(Viconst{k},2))*log(2*pi) + 0.5*log(abs(det(S0bar))) + ...
|
||||
0.5*log(abs(det(Spbar))) - 0.5*(bk'*S0bar*bk+(gk-gbark)'*Spbar*(gk-gbark));
|
||||
%--- For p(A+|Y,a0) only.
|
||||
tmpd = gk - Pmat{k}*bk;
|
||||
Apexpt = Apexpt - 0.5*tmpd'*(Hpinv{k}*tmpd);
|
||||
end
|
||||
vlog_a0p
|
||||
|
||||
%--------- The log value of p(Y|A0,A+) at some point such as the peak. ----------
|
||||
%--------- Note that logMarLHres is the same as vlog_Y_a, just to double check. ----------
|
||||
vlog_Y_a = -0.5*nvar*fss*log(2*pi) + fss*log(abs(det(A0xhat))) + Yexpt
|
||||
% a: given a0 and a+
|
||||
logMarLHres = 0; % Initialize log of the marginal likelihood (restricted or constant parameters).
|
||||
for ki=1:fss %ndobs+1:fss % Forward recursion to get the marginal likelihood. See F on p.19 and pp. 48-49.
|
||||
%---- Restricted log marginal likelihood function (constant parameters).
|
||||
[A0l,A0u] = lu(A0xhat);
|
||||
ada = sum(log(abs(diag(A0u)))); % log|A0|
|
||||
termexp = y(ki,:)*A0xhat - phi(ki,:)*Apxhat; % 1-by-nvar
|
||||
logMarLHres = logMarLHres - (0.5*nvar)*log(2*pi) + ada - 0.5*termexp*termexp'; % log value
|
||||
end
|
||||
logMarLHres
|
||||
|
||||
|
||||
%--------- The log value of p(A+|Y,A0) at some point such as the peak ----------
|
||||
totparsp = 0.0;
|
||||
tmpd = 0.0;
|
||||
for k=1:nvar
|
||||
totparsp = totparsp + size(Viconst{k},2);
|
||||
tmpd = tmpd + 0.5*log(abs(det(Hpinv{k})));
|
||||
end
|
||||
vlog_ap_Ya0 = -0.5*totparsp*log(2*pi) + tmpd + Apexpt;
|
||||
|
||||
|
||||
|
||||
|
||||
%===================================
|
||||
% Compute p(a0,k|Y,ao) at some point such as the peak (in this situation, we simply
|
||||
% generate results from the original Gibbs sampler). See FORECAST (2) pp.70-71
|
||||
%===================================
|
||||
%--- Global set up for Gibbs.
|
||||
[Tinv,UT] = fn_gibbsrvar_setup(H0inv, Uiconst, Hpinv, Pmat, Viconst, nvar, fss);
|
||||
%
|
||||
vlog_a0_Yao = zeros(nvar,1);
|
||||
% the log value of p(a0k|Y,ao) where ao: other a's at some point such as the peak of ONLY some a0's
|
||||
vlog=zeros(ndraws2,1);
|
||||
for k=1:nvar
|
||||
bk = Uiconst{k}'*A0xhat(:,k);
|
||||
indx_ks=[k:nvar]; % the columns that exclude 1-(k-1)th columns
|
||||
A0gbs0 = A0hat; % starting at some point such as the peak
|
||||
nk = n0(k);
|
||||
|
||||
if k<nvar
|
||||
%--------- The 1st set of draws to be tossed away. ------------------
|
||||
for draws = 1:ndraws1
|
||||
if ~mod(draws,nbuffer)
|
||||
disp(' ')
|
||||
disp(sprintf('The %dth column or equation in A0 with %d 1st tossed-away draws in Gibbs',k,draws))
|
||||
end
|
||||
A0gbs1 = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
|
||||
A0gbs0=A0gbs1; % repeat the Gibbs sampling
|
||||
end
|
||||
|
||||
|
||||
%--------- The 2nd set of draws to be used. ------------------
|
||||
for draws = 1:ndraws2
|
||||
if ~mod(draws,nbuffer)
|
||||
disp(' ')
|
||||
disp(sprintf('The %dth column or equation in A0 with %d usable draws in Gibbs',k,draws))
|
||||
end
|
||||
[A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
|
||||
%------ See p.71, Forecast (II).
|
||||
%------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
|
||||
Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II).
|
||||
gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
|
||||
[Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
|
||||
%
|
||||
vlog(draws) = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
|
||||
0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
|
||||
0.5*fss*bk'*(Vtr\(Vtr'\bk));
|
||||
|
||||
A0gbs0=A0gbs1; % repeat the Gibbs sampling
|
||||
end
|
||||
vlogm=max(vlog);
|
||||
qlog=vlog-vlogm;
|
||||
vlogxhat=vlogm-log(ndraws2)+log(sum(exp(qlog)));
|
||||
vlog_a0_Yao(k) = vlogxhat;
|
||||
% The log value of p(a0_k|Y,a_others) where a_others: other a's at some point such as the peak of ONLY some a0's
|
||||
else
|
||||
disp(' ')
|
||||
disp(sprintf('The last(6th) column or equation in A0 with no Gibbs draws'))
|
||||
[A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks)
|
||||
%------ See p.71, Forecast (II).
|
||||
%------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
|
||||
Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II).
|
||||
gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
|
||||
[Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
|
||||
%
|
||||
vloglast = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
|
||||
0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
|
||||
0.5*fss*bk'*(Vtr\(Vtr'\bk));
|
||||
vlog_a0_Yao(k) = vloglast;
|
||||
end
|
||||
end
|
||||
ndraws2
|
||||
|
||||
disp('Prior pdf -- log(p(a0hat, a+hat)):');
|
||||
vlog_a0p
|
||||
disp('LH pdf -- log(p(Y|a0hat, a+hat)):');
|
||||
vlog_Y_a
|
||||
disp('Posterior Kernal -- logp(ahat) + logp(Y|ahat):');
|
||||
vlog_Y_a + vlog_a0p
|
||||
disp('Posterior pdf -- log(p(a0_i_hat|a0_other_hat, Y)):');
|
||||
vlog_a0_Yao
|
||||
disp('Posterior pdf -- log(p(aphat|a0hat, Y)):');
|
||||
vlog_ap_Ya0
|
||||
|
||||
%--------- The value of marginal density p(Y) ----------
|
||||
disp(' ');
|
||||
disp(' ');
|
||||
disp('************ Marginal Likelihood of Y or Marginal Data Density: ************');
|
||||
vlogY = vlog_a0p+vlog_Y_a-sum(vlog_a0_Yao)-vlog_ap_Ya0
|
||||
|
|
|
|||
|
|
@ -1,162 +1,162 @@
|
|||
function ms_write_markov_file(fname, options)
|
||||
% function ms_write_markov_file(fname, options)
|
||||
%
|
||||
% INPUTS
|
||||
% fname: (string) name of markov file
|
||||
% options: (struct) options
|
||||
%
|
||||
% OUTPUTS
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
|
||||
% Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
n_chains = length(options.ms.ms_chain);
|
||||
nvars = size(options.varobs,1);
|
||||
|
||||
fh = fopen(fname,'w');
|
||||
%/******************************************************************************/
|
||||
%/********************* Markov State Variable Information **********************/
|
||||
%/******************************************************************************/
|
||||
|
||||
fprintf(fh,'//== Markov State Variables with Simple Restrictions ==//\n\n');
|
||||
|
||||
|
||||
%//This number is NOT used but read in.
|
||||
fprintf(fh,'//== Number Observations ==//\n');
|
||||
fprintf(fh,'0\n\n');
|
||||
|
||||
fprintf(fh,'//== Number independent state variables ==//\n');
|
||||
fprintf(fh,'%d\n\n',n_chains);
|
||||
|
||||
for i_chain = 1:n_chains
|
||||
|
||||
%//=====================================================//
|
||||
%//== state_variable[i] (1 <= i <= n_state_variables) ==//
|
||||
%//=====================================================//
|
||||
fprintf(fh,'//== Number of states for state_variable[%d] ==//\n', ...
|
||||
i_chain);
|
||||
n_states = length(options.ms.ms_chain(i_chain).regime);
|
||||
fprintf(fh,'%d\n\n',n_states);
|
||||
|
||||
%//== 03/15/06: DW TVBVAR code reads the data below and overwrite the prior data read somewhere else if any.
|
||||
%//== Each column contains the parameters for a Dirichlet prior on the corresponding
|
||||
%//== column of the transition matrix. Each element must be positive. For each column,
|
||||
%//== the relative size of the prior elements determine the relative size of the elements
|
||||
%//== of the transition matrix and overall larger sizes implies a tighter prior.
|
||||
fprintf(fh,['//== Transition matrix prior for state_variable[%d] ==//\n'], ...
|
||||
i_chain);
|
||||
Alpha = ones(n_states,n_states);
|
||||
for i_state = 1:n_states
|
||||
p = 1-1/options.ms.ms_chain(i_chain).regime(i_state).duration;
|
||||
Alpha(i_state,i_state) = p*(n_states-1)/(1-p);
|
||||
fprintf(fh,'%22.16f',Alpha(i_state,:));
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
|
||||
fprintf(fh,['\n//== Dirichlet dimensions for state_variable[%d] ' ...
|
||||
'==//\n'],i_chain);
|
||||
% fprintf(fh,'%d ',repmat(n_states,1,n_states));
|
||||
fprintf(fh,'%d ',repmat(2,1,n_states));
|
||||
fprintf(fh,'\n\n');
|
||||
|
||||
%//== The jth restriction matrix is n_states-by-free[j]. Each row of the restriction
|
||||
%//== matrix has exactly one non-zero entry and the sum of each column must be one.
|
||||
fprintf(fh,['//== Column restrictions for state_variable[%d] ' ...
|
||||
'==//\n'],i_chain);
|
||||
for i_state = 1:n_states
|
||||
if i_state == 1
|
||||
M = eye(n_states,2);
|
||||
elseif i_state == n_states
|
||||
M = [zeros(n_states-2,2); eye(2)];
|
||||
else
|
||||
M = zeros(n_states,2);
|
||||
M(i_state+[-1 1],1) = ones(2,1)/2;
|
||||
M(i_state,2) = 1;
|
||||
disp(M)
|
||||
end
|
||||
for j_state = 1:n_states
|
||||
fprintf(fh,'%f ',M(j_state,:));
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
end
|
||||
|
||||
%/******************************************************************************/
|
||||
%/******************************* VAR Parameters *******************************/
|
||||
%/******************************************************************************/
|
||||
%//NOT read
|
||||
fprintf(fh,'//== Number Variables ==//\n');
|
||||
fprintf(fh,'%d\n\n',nvars);
|
||||
|
||||
%//NOT read
|
||||
fprintf(fh,'//== Number Lags ==//\n');
|
||||
fprintf(fh,'%d\n\n',options.ms.nlags);
|
||||
|
||||
%//NOT read
|
||||
fprintf(fh,'//== Exogenous Variables ==//\n');
|
||||
fprintf(fh,'1\n\n');
|
||||
|
||||
|
||||
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
|
||||
%this state variable controls the jth column of A0 and Aplus
|
||||
fprintf(fh,['//== Controlling states variables for coefficients ==//\' ...
|
||||
'n']);
|
||||
|
||||
for i_var = 1:nvars
|
||||
for i_chain = 1:n_chains
|
||||
if ~isfield(options.ms.ms_chain(i_chain),'svar_coefficients') ...
|
||||
|| isempty(options.ms.ms_chain(i_chain).svar_coefficients)
|
||||
i_equations = 0;
|
||||
else
|
||||
i_equations = ...
|
||||
options.ms.ms_chain(i_chain).svar_coefficients.equations;
|
||||
end
|
||||
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
|
||||
fprintf(fh,'%d ',1);
|
||||
else
|
||||
fprintf(fh,'%d ',0);
|
||||
end
|
||||
end
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
|
||||
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
|
||||
%this state variable controls the jth diagonal element of Xi
|
||||
fprintf(fh,'\n//== Controlling states variables for variance ==//\n');
|
||||
for i_var = 1:nvars
|
||||
for i_chain = 1:n_chains
|
||||
if ~isfield(options.ms.ms_chain(i_chain),'svar_variances') ...
|
||||
|| isempty(options.ms.ms_chain(i_chain).svar_variances)
|
||||
i_equations = 0;
|
||||
else
|
||||
i_equations = ...
|
||||
options.ms.ms_chain(i_chain).svar_variances.equations;
|
||||
end
|
||||
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
|
||||
fprintf(fh,'%d ',1);
|
||||
else
|
||||
fprintf(fh,'%d ',0);
|
||||
end
|
||||
end
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
|
||||
fclose(fh);
|
||||
function ms_write_markov_file(fname, options)
|
||||
% function ms_write_markov_file(fname, options)
|
||||
%
|
||||
% INPUTS
|
||||
% fname: (string) name of markov file
|
||||
% options: (struct) options
|
||||
%
|
||||
% OUTPUTS
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
|
||||
% Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
n_chains = length(options.ms.ms_chain);
|
||||
nvars = size(options.varobs,1);
|
||||
|
||||
fh = fopen(fname,'w');
|
||||
%/******************************************************************************/
|
||||
%/********************* Markov State Variable Information **********************/
|
||||
%/******************************************************************************/
|
||||
|
||||
fprintf(fh,'//== Markov State Variables with Simple Restrictions ==//\n\n');
|
||||
|
||||
|
||||
%//This number is NOT used but read in.
|
||||
fprintf(fh,'//== Number Observations ==//\n');
|
||||
fprintf(fh,'0\n\n');
|
||||
|
||||
fprintf(fh,'//== Number independent state variables ==//\n');
|
||||
fprintf(fh,'%d\n\n',n_chains);
|
||||
|
||||
for i_chain = 1:n_chains
|
||||
|
||||
%//=====================================================//
|
||||
%//== state_variable[i] (1 <= i <= n_state_variables) ==//
|
||||
%//=====================================================//
|
||||
fprintf(fh,'//== Number of states for state_variable[%d] ==//\n', ...
|
||||
i_chain);
|
||||
n_states = length(options.ms.ms_chain(i_chain).regime);
|
||||
fprintf(fh,'%d\n\n',n_states);
|
||||
|
||||
%//== 03/15/06: DW TVBVAR code reads the data below and overwrite the prior data read somewhere else if any.
|
||||
%//== Each column contains the parameters for a Dirichlet prior on the corresponding
|
||||
%//== column of the transition matrix. Each element must be positive. For each column,
|
||||
%//== the relative size of the prior elements determine the relative size of the elements
|
||||
%//== of the transition matrix and overall larger sizes implies a tighter prior.
|
||||
fprintf(fh,['//== Transition matrix prior for state_variable[%d] ==//\n'], ...
|
||||
i_chain);
|
||||
Alpha = ones(n_states,n_states);
|
||||
for i_state = 1:n_states
|
||||
p = 1-1/options.ms.ms_chain(i_chain).regime(i_state).duration;
|
||||
Alpha(i_state,i_state) = p*(n_states-1)/(1-p);
|
||||
fprintf(fh,'%22.16f',Alpha(i_state,:));
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
|
||||
fprintf(fh,['\n//== Dirichlet dimensions for state_variable[%d] ' ...
|
||||
'==//\n'],i_chain);
|
||||
% fprintf(fh,'%d ',repmat(n_states,1,n_states));
|
||||
fprintf(fh,'%d ',repmat(2,1,n_states));
|
||||
fprintf(fh,'\n\n');
|
||||
|
||||
%//== The jth restriction matrix is n_states-by-free[j]. Each row of the restriction
|
||||
%//== matrix has exactly one non-zero entry and the sum of each column must be one.
|
||||
fprintf(fh,['//== Column restrictions for state_variable[%d] ' ...
|
||||
'==//\n'],i_chain);
|
||||
for i_state = 1:n_states
|
||||
if i_state == 1
|
||||
M = eye(n_states,2);
|
||||
elseif i_state == n_states
|
||||
M = [zeros(n_states-2,2); eye(2)];
|
||||
else
|
||||
M = zeros(n_states,2);
|
||||
M(i_state+[-1 1],1) = ones(2,1)/2;
|
||||
M(i_state,2) = 1;
|
||||
disp(M)
|
||||
end
|
||||
for j_state = 1:n_states
|
||||
fprintf(fh,'%f ',M(j_state,:));
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
end
|
||||
|
||||
%/******************************************************************************/
|
||||
%/******************************* VAR Parameters *******************************/
|
||||
%/******************************************************************************/
|
||||
%//NOT read
|
||||
fprintf(fh,'//== Number Variables ==//\n');
|
||||
fprintf(fh,'%d\n\n',nvars);
|
||||
|
||||
%//NOT read
|
||||
fprintf(fh,'//== Number Lags ==//\n');
|
||||
fprintf(fh,'%d\n\n',options.ms.nlags);
|
||||
|
||||
%//NOT read
|
||||
fprintf(fh,'//== Exogenous Variables ==//\n');
|
||||
fprintf(fh,'1\n\n');
|
||||
|
||||
|
||||
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
|
||||
%this state variable controls the jth column of A0 and Aplus
|
||||
fprintf(fh,['//== Controlling states variables for coefficients ==//\' ...
|
||||
'n']);
|
||||
|
||||
for i_var = 1:nvars
|
||||
for i_chain = 1:n_chains
|
||||
if ~isfield(options.ms.ms_chain(i_chain),'svar_coefficients') ...
|
||||
|| isempty(options.ms.ms_chain(i_chain).svar_coefficients)
|
||||
i_equations = 0;
|
||||
else
|
||||
i_equations = ...
|
||||
options.ms.ms_chain(i_chain).svar_coefficients.equations;
|
||||
end
|
||||
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
|
||||
fprintf(fh,'%d ',1);
|
||||
else
|
||||
fprintf(fh,'%d ',0);
|
||||
end
|
||||
end
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
|
||||
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
|
||||
%this state variable controls the jth diagonal element of Xi
|
||||
fprintf(fh,'\n//== Controlling states variables for variance ==//\n');
|
||||
for i_var = 1:nvars
|
||||
for i_chain = 1:n_chains
|
||||
if ~isfield(options.ms.ms_chain(i_chain),'svar_variances') ...
|
||||
|| isempty(options.ms.ms_chain(i_chain).svar_variances)
|
||||
i_equations = 0;
|
||||
else
|
||||
i_equations = ...
|
||||
options.ms.ms_chain(i_chain).svar_variances.equations;
|
||||
end
|
||||
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
|
||||
fprintf(fh,'%d ',1);
|
||||
else
|
||||
fprintf(fh,'%d ',0);
|
||||
end
|
||||
end
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
|
||||
fclose(fh);
|
||||
|
|
|
|||
|
|
@ -1,68 +1,68 @@
|
|||
function ms_write_mhm_input(fname, options_ms)
|
||||
% function ms_write_mhm_input(fname, options_ms)
|
||||
%
|
||||
% INPUTS
|
||||
% fname: (string) filename
|
||||
% options_ms: (struct) options
|
||||
%
|
||||
% OUTPUTS
|
||||
% none
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
|
||||
% Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
fh = fopen(fname,'w');
|
||||
|
||||
|
||||
fprintf(fh,'/**********************************************************\n');
|
||||
fprintf(fh,' *** This input file is read by swzmsbvar_mhm_1 and swzmsbvar_mhm_1.exe only, NOT by swzmsbvar_printdraws.exe.\n');
|
||||
fprintf(fh,' ***\n');
|
||||
fprintf(fh,' **********************************************************/\n');
|
||||
|
||||
fprintf(fh,'\n\n//------------- 1st set of posterior draws to find optimal scales for Metropolis (30000). ---------------\n');
|
||||
fprintf(fh,'//== number draws for first burn-in ==// //For determining the Metropolis scales only.\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_1);
|
||||
|
||||
fprintf(fh,'//------------- MCMC burn-in draws once the Metropolis scales (previous stage) are fixed. --------------\n');
|
||||
fprintf(fh,'//------------- 2nd set of standard burn-in posterior draws to throw away the initial draws (10000). ---------------\n');
|
||||
fprintf(fh,'//== number draws for second burn-in ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_2);
|
||||
|
||||
fprintf(fh,'//--------------- 1st set of posterior draws to compute the mean and variance for the weighting function in the MHM (200000) ----------------\n');
|
||||
fprintf(fh,'//== number draws to estimate mean and variance ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_mean_var_estimate);
|
||||
|
||||
fprintf(fh,'//--------------- Only applied to mhm_2 process: total number of MCMC draws = thinning factor * 2nd set of saved posterior draws ----------------\n');
|
||||
fprintf(fh,'//== thinning factor for modified harmonic mean process ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.thinning_factor);
|
||||
|
||||
fprintf(fh,'//--------------- 2nd set of saved posterior draws from MHM_2 (second stage): saved draws AFTER thinning (1000000) ----------------\n');
|
||||
fprintf(fh,'//== number draws for modified harmonic mean process ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_modified_harmonic_mean);
|
||||
|
||||
fprintf(fh,'//------- 1st stage: computing all three tightness factors for Dirichlet. ---------\n');
|
||||
fprintf(fh,'//------- 2nd stage: hard-code the second scale factor (in principle, we can do all three). ---------\n');
|
||||
fprintf(fh,'//------- It seems that Dan''s code only use the first element of the following scales. The scale applies to the Dirichlet''s hyperparameter alpha for the diagonal of the transition matrix in the weighting function. Note that the weighting function for the transition matrix parameters is Dirichlet. ---------\n');
|
||||
|
||||
fprintf(fh,'//== scale values for Dirichlet distribution ==//\n');
|
||||
fprintf(fh,'3\n\n');
|
||||
fprintf(fh,'%f ',options_ms.dirichlet_scale);
|
||||
fprintf(fh,'\n');
|
||||
function ms_write_mhm_input(fname, options_ms)
|
||||
% function ms_write_mhm_input(fname, options_ms)
|
||||
%
|
||||
% INPUTS
|
||||
% fname: (string) filename
|
||||
% options_ms: (struct) options
|
||||
%
|
||||
% OUTPUTS
|
||||
% none
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
|
||||
% Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
fh = fopen(fname,'w');
|
||||
|
||||
|
||||
fprintf(fh,'/**********************************************************\n');
|
||||
fprintf(fh,' *** This input file is read by swzmsbvar_mhm_1 and swzmsbvar_mhm_1.exe only, NOT by swzmsbvar_printdraws.exe.\n');
|
||||
fprintf(fh,' ***\n');
|
||||
fprintf(fh,' **********************************************************/\n');
|
||||
|
||||
fprintf(fh,'\n\n//------------- 1st set of posterior draws to find optimal scales for Metropolis (30000). ---------------\n');
|
||||
fprintf(fh,'//== number draws for first burn-in ==// //For determining the Metropolis scales only.\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_1);
|
||||
|
||||
fprintf(fh,'//------------- MCMC burn-in draws once the Metropolis scales (previous stage) are fixed. --------------\n');
|
||||
fprintf(fh,'//------------- 2nd set of standard burn-in posterior draws to throw away the initial draws (10000). ---------------\n');
|
||||
fprintf(fh,'//== number draws for second burn-in ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_2);
|
||||
|
||||
fprintf(fh,'//--------------- 1st set of posterior draws to compute the mean and variance for the weighting function in the MHM (200000) ----------------\n');
|
||||
fprintf(fh,'//== number draws to estimate mean and variance ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_mean_var_estimate);
|
||||
|
||||
fprintf(fh,'//--------------- Only applied to mhm_2 process: total number of MCMC draws = thinning factor * 2nd set of saved posterior draws ----------------\n');
|
||||
fprintf(fh,'//== thinning factor for modified harmonic mean process ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.thinning_factor);
|
||||
|
||||
fprintf(fh,'//--------------- 2nd set of saved posterior draws from MHM_2 (second stage): saved draws AFTER thinning (1000000) ----------------\n');
|
||||
fprintf(fh,'//== number draws for modified harmonic mean process ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_modified_harmonic_mean);
|
||||
|
||||
fprintf(fh,'//------- 1st stage: computing all three tightness factors for Dirichlet. ---------\n');
|
||||
fprintf(fh,'//------- 2nd stage: hard-code the second scale factor (in principle, we can do all three). ---------\n');
|
||||
fprintf(fh,'//------- It seems that Dan''s code only use the first element of the following scales. The scale applies to the Dirichlet''s hyperparameter alpha for the diagonal of the transition matrix in the weighting function. Note that the weighting function for the transition matrix parameters is Dirichlet. ---------\n');
|
||||
|
||||
fprintf(fh,'//== scale values for Dirichlet distribution ==//\n');
|
||||
fprintf(fh,'3\n\n');
|
||||
fprintf(fh,'%f ',options_ms.dirichlet_scale);
|
||||
fprintf(fh,'\n');
|
||||
fclose(fh);
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,166 +1,166 @@
|
|||
%function []= msstart_setup(options_)
|
||||
|
||||
% Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
% ** ONLY UNDER UNIX SYSTEM
|
||||
%path(path,'/usr2/f1taz14/mymatlab')
|
||||
|
||||
|
||||
|
||||
%===========================================
|
||||
% Exordium I
|
||||
%===========================================
|
||||
format short g % format
|
||||
%
|
||||
%options_.ms.freq = 4; % quarters or months
|
||||
%options_.ms.initial_year=1959; % beginning of the year
|
||||
%options_.ms.initial_subperiod=1; % begining of the quarter or month
|
||||
%options_.ms.final_year=2005; % final year
|
||||
%options_.ms.final_subperiod=4; % final month or quarter
|
||||
nData=(options_.ms.final_year-options_.ms.initial_year)*options_.ms.freq + (options_.ms.final_subperiod-options_.ms.initial_subperiod+1);
|
||||
% total number of the available data -- this is all you have
|
||||
|
||||
%*** Load data and series
|
||||
%load datainf_argen.prn % the default name for the variable is "options_.ms.data".
|
||||
%load datacbogdpffr.prn
|
||||
%options_.ms.data = datacbogdpffr;
|
||||
%clear datacbogdpffr;
|
||||
[nt,ndv]=size(options_.data);
|
||||
if nt~=nData
|
||||
disp(' ')
|
||||
warning(sprintf('nt=%d, Caution: not equal to the length in the data',nt));
|
||||
%disp(sprintf('nt=%d, Caution: not equal to the length in the data',nt));
|
||||
disp('Press ctrl-c to abort')
|
||||
return
|
||||
end
|
||||
%--------
|
||||
%1 CBO output gap -- log(x_t)-log(x_t potential)
|
||||
%2 GDP deflator -- (P_t/P_{t-1})^4-1.0
|
||||
%2 FFR/100.
|
||||
options_.ms.vlist = [1:size(options_.varobs,1)]; % 1: U; 4: PCE inflation.
|
||||
options_.ms.varlist=cellstr(options_.varobs);
|
||||
options_.ms.log_var = sort(varlist_indices(options_.ms.vlistlog,options_.varobs)); % subset of "options_.ms.vlist. Variables in log level so that differences are in **monthly** growth, unlike R and U which are in annual percent (divided by 100 already).
|
||||
options_.ms.percent_var =setdiff(options_.ms.vlist,options_.ms.log_var);
|
||||
%options_.ms.restriction_fname='ftd_upperchol3v'; %Only used by msstart2.m.
|
||||
ylab = options_.ms.varlist;
|
||||
xlab = options_.ms.varlist;
|
||||
|
||||
%----------------
|
||||
nvar = size(options_.varobs,1); % number of endogenous variables
|
||||
nlogeno = length(options_.ms.log_var); % number of endogenous variables in options_.ms.log_var
|
||||
npereno = length(options_.ms.percent_var); % number of endogenous variables in options_.ms.percent_var
|
||||
if (nvar~=(nlogeno+npereno))
|
||||
disp(' ')
|
||||
warning('Check xlab, nlogeno or npereno to make sure of endogenous variables in options_.ms.vlist')
|
||||
disp('Press ctrl-c to abort')
|
||||
return
|
||||
elseif (nvar==length(options_.ms.vlist))
|
||||
nexo=1; % only constants as an exogenous variable. The default setting.
|
||||
elseif (nvar<length(options_.ms.vlist))
|
||||
nexo=length(options_.ms.vlist)-nvar+1;
|
||||
else
|
||||
disp(' ')
|
||||
warning('Make sure there are only nvar endogenous variables in options_.ms.vlist')
|
||||
disp('Press ctrl-c to abort')
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
%------- A specific sample is considered for estimation -------
|
||||
yrStart=options_.ms.initial_year;
|
||||
qmStart=options_.ms.initial_subperiod;
|
||||
yrEnd=options_.ms.final_year;
|
||||
qmEnd=options_.ms.final_subperiod;
|
||||
%options_.forecast = 4; % number of years for forecasting
|
||||
if options_.forecast<1
|
||||
error('To be safe, the number of forecast years should be at least 1')
|
||||
end
|
||||
ystr=num2str(yrEnd);
|
||||
forelabel = [ ystr(3:4) ':' num2str(qmEnd) ' Forecast'];
|
||||
|
||||
nSample=(yrEnd-yrStart)*options_.ms.freq + (qmEnd-qmStart+1);
|
||||
if qmEnd==options_.ms.freq
|
||||
E1yrqm = [yrEnd+1 1]; % first year and quarter (month) after the sample
|
||||
else
|
||||
E1yrqm = [yrEnd qmEnd+1]; % first year and quarter (month) after the sample
|
||||
end
|
||||
E2yrqm = [yrEnd+options_.forecast qmEnd]; % end at the last month (quarter) of a calendar year after the sample
|
||||
[fdates,nfqm]=fn_calyrqm(options_.ms.freq,E1yrqm,E2yrqm); % forecast dates and number of forecast dates
|
||||
[sdates,nsqm] = fn_calyrqm(options_.ms.freq,[yrStart qmStart],[yrEnd qmEnd]);
|
||||
% sdates: dates for the whole sample (including options_.ms.nlags)
|
||||
if nSample~=nsqm
|
||||
warning('Make sure that nSample is consistent with the size of sdates')
|
||||
disp('Hit any key to continue, or ctrl-c to abort')
|
||||
pause
|
||||
end
|
||||
imstp = 4*options_.ms.freq; % <<>> impulse responses (4 years)
|
||||
nayr = 4; %options_.forecast; % number of years before forecasting for plotting.
|
||||
|
||||
|
||||
%------- Prior, etc. -------
|
||||
%options_.ms.nlags = 4; % number of options_.ms.nlags
|
||||
%options_.ms.cross_restrictions = 0; % 1: cross-A0-and-A+ restrictions; 0: options_.ms.restriction_fname is all we have
|
||||
% Example for indxOres==1: restrictions of the form P(t) = P(t-1).
|
||||
%options_.ms.contemp_reduced_form = 0; % 1: contemporaneous recursive reduced form; 0: restricted (non-recursive) form
|
||||
%options_.ms.real_pseudo_forecast = 0; % 1: options_.ms.real_pseudo_forecast forecasts; 0: real time forecasts
|
||||
%options_.ms.bayesian_prior = 1; % 1: Bayesian prior; 0: no prior
|
||||
indxDummy = options_.ms.bayesian_prior; % 1: add dummy observations to the data; 0: no dummy added.
|
||||
%options_.ms.dummy_obs = 0; % No dummy observations for xtx, phi, fss, xdatae, etc. Dummy observations are used as an explicit prior in fn_rnrprior_covres_dobs.m.
|
||||
%if indxDummy
|
||||
% options_.ms.dummy_obs=nvar+1; % number of dummy observations
|
||||
%else
|
||||
% options_.ms.dummy_obs=0; % no dummy observations
|
||||
%end
|
||||
%=== The following mu is effective only if options_.ms.bayesian_prior==1.
|
||||
|
||||
mu = options_.ms.coefficients_prior_hyperparameters;
|
||||
|
||||
% mu(1): overall tightness and also for A0;
|
||||
% mu(2): relative tightness for A+;
|
||||
% mu(3): relative tightness for the constant term;
|
||||
% mu(4): tightness on lag decay; (1)
|
||||
% mu(5): weight on nvar sums of coeffs dummy observations (unit roots);
|
||||
% mu(6): weight on single dummy initial observation including constant
|
||||
% (cointegration, unit roots, and stationarity);
|
||||
%
|
||||
%
|
||||
hpmsmd = [0.0; 0.0];
|
||||
indxmsmdeqn = [0; 0; 0; 0]; %This option disenable using this in fn_rnrprior_covres_dobs.m
|
||||
|
||||
|
||||
tdf = 3; % degrees of freedom for t-dist for initial draw of the MC loop
|
||||
nbuffer = 1000; % a block or buffer of draws (buffer) that is saved to the disk (not memory)
|
||||
ndraws1=1*nbuffer; % 1st part of Monte Carlo draws
|
||||
ndraws2=10*ndraws1; % 2nd part of Monte Carlo draws
|
||||
% seednumber = options_.DynareRandomStreams.seed; %7910; %472534; % if 0, random state at each clock time
|
||||
% % good one 420 for [29 45], [29 54]
|
||||
% if seednumber
|
||||
% randn('state',seednumber);
|
||||
% rand('state',seednumber);
|
||||
% else
|
||||
% randn('state',fix(100*sum(clock)));
|
||||
% rand('state',fix(100*sum(clock)));
|
||||
% end
|
||||
% nstarts=1 % number of starting points
|
||||
% imndraws = nstarts*ndraws2; % total draws for impulse responses or forecasts
|
||||
%<<<<<<<<<<<<<<<<<<<
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
%function []= msstart_setup(options_)
|
||||
|
||||
% Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
% ** ONLY UNDER UNIX SYSTEM
|
||||
%path(path,'/usr2/f1taz14/mymatlab')
|
||||
|
||||
|
||||
|
||||
%===========================================
|
||||
% Exordium I
|
||||
%===========================================
|
||||
format short g % format
|
||||
%
|
||||
%options_.ms.freq = 4; % quarters or months
|
||||
%options_.ms.initial_year=1959; % beginning of the year
|
||||
%options_.ms.initial_subperiod=1; % begining of the quarter or month
|
||||
%options_.ms.final_year=2005; % final year
|
||||
%options_.ms.final_subperiod=4; % final month or quarter
|
||||
nData=(options_.ms.final_year-options_.ms.initial_year)*options_.ms.freq + (options_.ms.final_subperiod-options_.ms.initial_subperiod+1);
|
||||
% total number of the available data -- this is all you have
|
||||
|
||||
%*** Load data and series
|
||||
%load datainf_argen.prn % the default name for the variable is "options_.ms.data".
|
||||
%load datacbogdpffr.prn
|
||||
%options_.ms.data = datacbogdpffr;
|
||||
%clear datacbogdpffr;
|
||||
[nt,ndv]=size(options_.data);
|
||||
if nt~=nData
|
||||
disp(' ')
|
||||
warning(sprintf('nt=%d, Caution: not equal to the length in the data',nt));
|
||||
%disp(sprintf('nt=%d, Caution: not equal to the length in the data',nt));
|
||||
disp('Press ctrl-c to abort')
|
||||
return
|
||||
end
|
||||
%--------
|
||||
%1 CBO output gap -- log(x_t)-log(x_t potential)
|
||||
%2 GDP deflator -- (P_t/P_{t-1})^4-1.0
|
||||
%2 FFR/100.
|
||||
options_.ms.vlist = [1:size(options_.varobs,1)]; % 1: U; 4: PCE inflation.
|
||||
options_.ms.varlist=cellstr(options_.varobs);
|
||||
options_.ms.log_var = sort(varlist_indices(options_.ms.vlistlog,options_.varobs)); % subset of "options_.ms.vlist. Variables in log level so that differences are in **monthly** growth, unlike R and U which are in annual percent (divided by 100 already).
|
||||
options_.ms.percent_var =setdiff(options_.ms.vlist,options_.ms.log_var);
|
||||
%options_.ms.restriction_fname='ftd_upperchol3v'; %Only used by msstart2.m.
|
||||
ylab = options_.ms.varlist;
|
||||
xlab = options_.ms.varlist;
|
||||
|
||||
%----------------
|
||||
nvar = size(options_.varobs,1); % number of endogenous variables
|
||||
nlogeno = length(options_.ms.log_var); % number of endogenous variables in options_.ms.log_var
|
||||
npereno = length(options_.ms.percent_var); % number of endogenous variables in options_.ms.percent_var
|
||||
if (nvar~=(nlogeno+npereno))
|
||||
disp(' ')
|
||||
warning('Check xlab, nlogeno or npereno to make sure of endogenous variables in options_.ms.vlist')
|
||||
disp('Press ctrl-c to abort')
|
||||
return
|
||||
elseif (nvar==length(options_.ms.vlist))
|
||||
nexo=1; % only constants as an exogenous variable. The default setting.
|
||||
elseif (nvar<length(options_.ms.vlist))
|
||||
nexo=length(options_.ms.vlist)-nvar+1;
|
||||
else
|
||||
disp(' ')
|
||||
warning('Make sure there are only nvar endogenous variables in options_.ms.vlist')
|
||||
disp('Press ctrl-c to abort')
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
%------- A specific sample is considered for estimation -------
|
||||
yrStart=options_.ms.initial_year;
|
||||
qmStart=options_.ms.initial_subperiod;
|
||||
yrEnd=options_.ms.final_year;
|
||||
qmEnd=options_.ms.final_subperiod;
|
||||
%options_.forecast = 4; % number of years for forecasting
|
||||
if options_.forecast<1
|
||||
error('To be safe, the number of forecast years should be at least 1')
|
||||
end
|
||||
ystr=num2str(yrEnd);
|
||||
forelabel = [ ystr(3:4) ':' num2str(qmEnd) ' Forecast'];
|
||||
|
||||
nSample=(yrEnd-yrStart)*options_.ms.freq + (qmEnd-qmStart+1);
|
||||
if qmEnd==options_.ms.freq
|
||||
E1yrqm = [yrEnd+1 1]; % first year and quarter (month) after the sample
|
||||
else
|
||||
E1yrqm = [yrEnd qmEnd+1]; % first year and quarter (month) after the sample
|
||||
end
|
||||
E2yrqm = [yrEnd+options_.forecast qmEnd]; % end at the last month (quarter) of a calendar year after the sample
|
||||
[fdates,nfqm]=fn_calyrqm(options_.ms.freq,E1yrqm,E2yrqm); % forecast dates and number of forecast dates
|
||||
[sdates,nsqm] = fn_calyrqm(options_.ms.freq,[yrStart qmStart],[yrEnd qmEnd]);
|
||||
% sdates: dates for the whole sample (including options_.ms.nlags)
|
||||
if nSample~=nsqm
|
||||
warning('Make sure that nSample is consistent with the size of sdates')
|
||||
disp('Hit any key to continue, or ctrl-c to abort')
|
||||
pause
|
||||
end
|
||||
imstp = 4*options_.ms.freq; % <<>> impulse responses (4 years)
|
||||
nayr = 4; %options_.forecast; % number of years before forecasting for plotting.
|
||||
|
||||
|
||||
%------- Prior, etc. -------
|
||||
%options_.ms.nlags = 4; % number of options_.ms.nlags
|
||||
%options_.ms.cross_restrictions = 0; % 1: cross-A0-and-A+ restrictions; 0: options_.ms.restriction_fname is all we have
|
||||
% Example for indxOres==1: restrictions of the form P(t) = P(t-1).
|
||||
%options_.ms.contemp_reduced_form = 0; % 1: contemporaneous recursive reduced form; 0: restricted (non-recursive) form
|
||||
%options_.ms.real_pseudo_forecast = 0; % 1: options_.ms.real_pseudo_forecast forecasts; 0: real time forecasts
|
||||
%options_.ms.bayesian_prior = 1; % 1: Bayesian prior; 0: no prior
|
||||
indxDummy = options_.ms.bayesian_prior; % 1: add dummy observations to the data; 0: no dummy added.
|
||||
%options_.ms.dummy_obs = 0; % No dummy observations for xtx, phi, fss, xdatae, etc. Dummy observations are used as an explicit prior in fn_rnrprior_covres_dobs.m.
|
||||
%if indxDummy
|
||||
% options_.ms.dummy_obs=nvar+1; % number of dummy observations
|
||||
%else
|
||||
% options_.ms.dummy_obs=0; % no dummy observations
|
||||
%end
|
||||
%=== The following mu is effective only if options_.ms.bayesian_prior==1.
|
||||
|
||||
mu = options_.ms.coefficients_prior_hyperparameters;
|
||||
|
||||
% mu(1): overall tightness and also for A0;
|
||||
% mu(2): relative tightness for A+;
|
||||
% mu(3): relative tightness for the constant term;
|
||||
% mu(4): tightness on lag decay; (1)
|
||||
% mu(5): weight on nvar sums of coeffs dummy observations (unit roots);
|
||||
% mu(6): weight on single dummy initial observation including constant
|
||||
% (cointegration, unit roots, and stationarity);
|
||||
%
|
||||
%
|
||||
hpmsmd = [0.0; 0.0];
|
||||
indxmsmdeqn = [0; 0; 0; 0]; %This option disenable using this in fn_rnrprior_covres_dobs.m
|
||||
|
||||
|
||||
tdf = 3; % degrees of freedom for t-dist for initial draw of the MC loop
|
||||
nbuffer = 1000; % a block or buffer of draws (buffer) that is saved to the disk (not memory)
|
||||
ndraws1=1*nbuffer; % 1st part of Monte Carlo draws
|
||||
ndraws2=10*ndraws1; % 2nd part of Monte Carlo draws
|
||||
% seednumber = options_.DynareRandomStreams.seed; %7910; %472534; % if 0, random state at each clock time
|
||||
% % good one 420 for [29 45], [29 54]
|
||||
% if seednumber
|
||||
% randn('state',seednumber);
|
||||
% rand('state',seednumber);
|
||||
% else
|
||||
% randn('state',fix(100*sum(clock)));
|
||||
% rand('state',fix(100*sum(clock)));
|
||||
% end
|
||||
% nstarts=1 % number of starting points
|
||||
% imndraws = nstarts*ndraws2; % total draws for impulse responses or forecasts
|
||||
%<<<<<<<<<<<<<<<<<<<
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue