diff --git a/matlab/ms-sbvar/ms_mardd.m b/matlab/ms-sbvar/ms_mardd.m
index 66e11222b..078bd5b92 100644
--- a/matlab/ms-sbvar/ms_mardd.m
+++ b/matlab/ms-sbvar/ms_mardd.m
@@ -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 .
-
-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.
+
+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.
-
- 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 .
+
+ 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);
diff --git a/matlab/ms-sbvar/ms_write_mhm_input.m b/matlab/ms-sbvar/ms_write_mhm_input.m
index b5ec084a1..f9f560135 100644
--- a/matlab/ms-sbvar/ms_write_mhm_input.m
+++ b/matlab/ms-sbvar/ms_write_mhm_input.m
@@ -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 .
-
- 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 .
+
+ 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);
\ No newline at end of file
diff --git a/matlab/ms-sbvar/msstart2.m b/matlab/ms-sbvar/msstart2.m
index 48990a51c..14bb31554 100644
--- a/matlab/ms-sbvar/msstart2.m
+++ b/matlab/ms-sbvar/msstart2.m
@@ -1,808 +1,808 @@
-%function msstart2(options_)
-% This .m file is called for point graphics or error bands and
-% It starts for both data and Bayesian estimation (when IxEstimat==0,
-% no estimation but only data analysis), which allows you to set
-% (a) available data range,
-% (b) sample range,
-% (c) rearrangement of actual data such as mlog, qg, yg
-% (d) conditions of shocks 'Cms',
-% (c) conditions of variables 'nconstr',
-% (e) soft conditions 'nbancon,'
-% (f) produce point conditional forecast (at least conditions on variables).
-%
-% February 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 .
-
-% ** ONLY UNDER UNIX SYSTEM
-%path(path,'/usr2/f1taz14/mymatlab')
-
-%addpath('C:\SoftWDisk\MATLAB6p5\toolbox\cstz')
-
-
-msstart_setup
-%===========================================
-% Exordium II
-%===========================================
-%options_.ms.ncsk = 0; % conditional directly on shoocks. Unlike Cms, not on variables that back out shocks
-%options_.ms.nstd = 6; % number of standard deviations to cover the range in which distributions are put to bins
-%options_.ms.ninv = 1000; % the number of bins for grouping, say, impulse responses
-%options_.ms.Indxcol = [1:nvar]; % a vector of random columns in which MC draws are made.
-%
-%options_.ms.indxparr = 1; % 1: parameters random; 0: no randomness in parameters
- % Note, when 0, there is no effect from the values of options_.ms.IndxAp, options_.ms.Aband, etc.
-%options_.ms.indxovr = 0; % 1: distributions for other variables of interest; 0: no distribution.
- % Example: joint distribution of a(1) and a(2). Only for specific purposes
-%options_.ms.Aband = 1; % 1: error bands with only A0 and A+ random.
-%options_.ms.IndxAp = 1; % 1: generate draws of A+; 0: no such draws.
- % Note: when options_.ms.IndxAp=0, there is no effect from the values of options_.ms.options_.ms.options_.ms.options_.ms.indximf, IndxFore,
- % or options_.ms.apband.
-%options_.ms.apband = 1; % 1: error bands for A+; 0: no error bands for A+.
-%*** The following (impulse responses and forecasts) is used only if options_.ms.IndxAp=1
-%options_.ms.indximf = 1; % 1: generate draws of impulse responses; 0: no such draws (thus no effect
- % from options_.ms.imfband)
-%options_.ms.imfband = 1; % 1: error bands for impulse responses; 0: no error bands
-%options_.ms.indxfore = 0; % 1: generate draws of forecasts; 0: no such draws (thus no effect from options_.ms.foreband)
-%options_.ms.foreband = 0; % 1: error bands for out-of-sample forecasts; 0: no error bands
-%
-%options_.ms.indxgforhat = 1; % 1: plot unconditoinal forecasts; 0: no such plot
-rnum = nvar; % number of rows in the graph
-cnum = 1; % number of columns in the graph
-if rnum*cnum>
-nconstr=0 ; %6*nconstr1;
-options_.ms.eq_ms = []; % location of MS equation; if [], all shocks
-PorR = [4*ones(nconstr1,1);2*ones(nconstr1,1);3*ones(nconstr1,1)]; % the variable conditioned. 1: Pcm; 3: FFR; 4: CPI
-PorR = [PorR;1*ones(nconstr1,1);5*ones(nconstr1,1);6*ones(nconstr1,1)];
-%PorR = [3 5]; % the variable conditioned. 3: FFR; 5: CPI
-%PorR = 3;
-
-%%%%----------------------------------------
-% Conditions directly on future shocks
-%
-%options_.ms.cms = 0 % 1: condition on ms shocks; 0: disable this and "fidcnderr.m" gives
- % unconditional forecasts if nconstr = 0 as well; <<>>
-%options_.ms.ncms = 0; % number of the stance of policy; 0 if no tightening or loosening
-%options_.ms.eq_cms = 1; % location of MS shocks
-options_.ms.tlindx = 1*ones(1,options_.ms.ncms); % 1-by-options_.ms.ncms vector; 1: tightening; 0: loosen
-options_.ms.tlnumber = [0.5 0.5 0 0]; %94:4 % [2 2 1.5 1.5]; %79:9 %[1.5 1.5 1 1]; 90:9
- % 1-by-options_.ms.ncms vector; cut-off point for MS shocks
-TLmean = zeros(1,options_.ms.ncms);
- % unconditional, i.e., 0 mean, for the final report in the paper
-if options_.ms.cms
- options_.ms.eq_ms = [];
- % At least at this point, it makes no sense to have DLS type of options_.ms.eq_ms; 10/12/98
- if all(isfinite(options_.ms.tlnumber))
- for k=1:options_.ms.ncms
- TLmean(k) = lcnmean(options_.ms.tlnumber(k),options_.ms.tlindx(k));
- % shock mean magnitude. 1: tight; 0: loose
- % Never used for any subsequent computation but
- % simply used for the final report in the paper.
- %options_.ms.tlnumber(k) = fzero('lcutoff',0,[],[],TLmean(k))
- % get an idea about the cutoff point given TLmean instead
-
- end
- end
-else
- options_.ms.ncms = 0; % only for the use of the graph by msprobg.m
- options_.ms.tlnumber = NaN*ones(1,options_.ms.ncms);
- % -infinity, only for the use of the graph by msprobg.m
-end
-
-
-%%%%----------------------------------------
-% Soft conditions on variables
-%
-%cnum = 0 % # of band condtions; when 0, disable this option
- % Note (different from "fidencon") that each condition corres. to variable
-%options_.ms.banact = 1; % 1: use infor on actual; 0: preset without infor on actual
-if cnum
- banindx = cell(cnum,1); % index for each variable or conditon
- banstp = cell(cnum,1); % steps: annual in general
- banvar = zeros(cnum,1); % varables: annual in general
- banval = cell(cnum,1); % band value (each variable occupy a cell)
- badval{1} = zeros(length(banstp{1}),2); % 2: lower or higher bound
-
- banstp{1} = 1:4; % 3 or 4 years
- banvar(1) = 3; % 3: FFR; 5: CPI
- if ~options_.ms.banact
- for i=1:length(banstp{1})
- banval{1}(i,:) = [5.0 10.0];
- end
- end
-end
-%
-pause(1)
-disp(' ')
-disp('For uncondtional forecasts, set nconstr = options_.ms.cms = cnum = 0')
-pause(1)
-%
-%=================================================
-%====== End of exordium ==========================
-%=================================================
-
-
-
-
-
-%(1)--------------------------------------
-% Further data analysis
-%(1)--------------------------------------
-%
-if (options_.ms.freq==12)
- nStart=(yrStart-options_.ms.initial_year )*12+qmStart-options_.ms.initial_subperiod ; % positive number of months at the start
- nEnd=(yrEnd-options_.ms.final_year )*12+qmEnd-options_.ms.final_subperiod ; % negative number of months towards end
-elseif (options_.ms.freq==4)
- nStart=(yrStart-options_.ms.initial_year )*4+qmStart-options_.ms.initial_subperiod ; % positive number of months at the start
- nEnd=(yrEnd-options_.ms.final_year )*4+qmEnd-options_.ms.final_subperiod ; % negative number of months towards end
-else
- error('Error: this code is only good for monthly/quarterly data!!!')
- return
-end
-%
-if nEnd>0 || nStart<0
- disp('Warning: this particular sample consider is out of bounds of the data!!!')
- return
-end
-%*** Note, both xdgel and xdata have the same start with the specific sample
-xdgel=options_.data(nStart+1:nData+nEnd,options_.ms.vlist);
- % gel: general options_.data within sample (nSample)
-if ~(nSample==size(xdgel,1))
- warning('The sample size (including options_.ms.nlags ) and data are incompatible')
- disp('Check to make sure nSample and size(xdgel,1) are the same')
- return
-end
-%
-baddata = find(isnan(xdgel));
-if ~isempty(baddata)
- warning('Some data for this selected sample are actually unavailable.')
- disp('Hit any key to continue, or ctrl-c to abort')
- pause
-end
-%
-if options_.ms.initial_subperiod ==1
- yrB = options_.ms.initial_year ; qmB = options_.ms.initial_subperiod ;
-else
- yrB = options_.ms.initial_year +1; qmB = 1;
-end
-yrF = options_.ms.final_year ; qmF = options_.ms.final_subperiod ;
-[Mdate,tmp] = fn_calyrqm(options_.ms.freq,[options_.ms.initial_year options_.ms.initial_subperiod ],[options_.ms.final_year options_.ms.final_subperiod ]);
-xdatae=[Mdate options_.data(1:nData,options_.ms.vlist)];
- % beyond sample into forecast horizon until the end of the data options_.ms.final_year :options_.ms.final_subperiod
- % Note: may contain NaN data. So must be careful about its use
-
-%=========== Obtain prior-period, period-to-last period, and annual growth rates
-[yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(xdatae,options_.ms.freq,options_.ms.log_var,options_.ms.percent_var,[yrB qmB],[yrF qmF]);
-qdates = zeros(size(yactqmyge,1),1);
-for ki=1:length(qdates)
- qdates(ki) = yactqmyge(1,1) + (yactqmyge(1,2)+ki-2)/options_.ms.freq;
-end
-for ki=1:nvar
- figure
- plot(qdates, yactqmyge(:,2+ki)/100)
- xlabel(options_.ms.varlist{ki})
-end
-save outactqmygdata.prn yactqmyge -ascii
-
-
-
-%=========== Write the output on the screen or to a file in an organized way ==============
-%disp([sprintf('%4.0f %2.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',yactyrge')])
-spstr1 = 'disp([sprintf(';
-spstr2 = '%4.0f %2.0f';
-yactyrget=yactyrge';
-for ki=1:length(options_.ms.vlist)
- if ki==length(options_.ms.vlist)
- spstr2 = [spstr2 ' %8.3f\n'];
- else
- spstr2 = [spstr2 ' %8.3f'];
- end
-end
-spstr = [spstr1 'spstr2' ', yactyrget)])'];
-eval(spstr)
-
-%
-fid = fopen('outyrqm.prn','w');
-%fprintf(fid,'%4.0f %2.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',yactyrge');
-fpstr1 = 'fprintf(fid,';
-fpstr2 = '%4.0f %2.0f';
-for ki=1:nvar
- if ki==nvar
- fpstr2 = [fpstr2 ' %8.3f\n'];
- else
- fpstr2 = [fpstr2 ' %8.3f'];
- end
-end
-fpstr = [fpstr1 'fpstr2' ', yactyrget);'];
-eval(fpstr)
-fclose(fid);
-
-
-
-if options_.ms.indxestima
- %(2)----------------------------------------------------------------------------
- % Estimation
- % ML forecast and impulse responses
- % Hard or soft conditions for conditional forecasts
- %(2)----------------------------------------------------------------------------
- %
- %* Arranged data information, WITHOUT dummy obs when 0 after mu is used. See fn_rnrprior_covres_dobs.m for using the dummy
- % observations as part of an explicit prior.
- [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh] = fn_dataxy(nvar,options_.ms.nlags ,xdgel,mu,0,nexo);
- if qmStart+options_.ms.nlags -options_.ms.dummy_obs >0
- qmStartEsti = rem(qmStart+options_.ms.nlags -options_.ms.dummy_obs ,options_.ms.freq); % dummy observations are included in the sample.
- if (~qmStartEsti)
- qmStartEsti = options_.ms.freq;
- end
- yrStartEsti = yrStart + floor((qmStart+options_.ms.nlags -options_.ms.dummy_obs )/(options_.ms.freq+0.01));
- % + 0.01 (or any number < 1) is used so that qmStart+options_.ms.nlags -options_.ms.dummy_obs ==?*options_.ms.freq doesn't give us an extra year forward.
- else
- qmStartEsti = options_.ms.freq + rem(qmStart+options_.ms.nlags -options_.ms.dummy_obs ,options_.ms.freq); % dummy observations are included in the sample.
- if (qmStart+options_.ms.nlags -options_.ms.dummy_obs ==0)
- yrStartEsti = yrStart - 1; % one year back.
- else
- yrStartEsti = yrStart + floor((qmStart+options_.ms.nlags -options_.ms.dummy_obs )/(options_.ms.freq-0.01));
- % - 0.01 (or any number < 1) is used so that qmStart+options_.ms.nlags -options_.ms.dummy_obs ==-?*options_.ms.freq give us an extra year back.
- end
- end
- dateswd = fn_dataext([yrStartEsti qmStartEsti],[yrEnd qmEnd],xdatae(:,[1:2])); % dates with dummies
- phie = [dateswd phi];
- ye = [dateswd y];
-
- %* Obtain linear restrictions
- [Uiconst,Viconst,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,nvar,nexo,options_.ms );
- if min(n0)==0
- disp(' ')
- warning('A0: restrictions in dlrprior.m give no free parameter in one of equations')
- disp('Press ctrl-c to abort')
- pause
- elseif min(np)==0
- disp(' ')
- warning('Ap: Restrictions in dlrprior.m give no free parameter in one of equations')
- disp('Press ctrl-c to abort')
- pause
- end
-
- if options_.ms.contemp_reduced_form
- Uiconst=cell(nvar,1); Viconst=cell(ncoef,1);
- for kj=1:nvar
- Uiconst{kj} = eye(nvar); Viconst{kj} = eye(ncoef);
- end
- end
-
- if options_.ms.bayesian_prior
- %*** Obtains asymmetric prior (with no linear restrictions) with dummy observations as part of an explicit prior (i.e,
- % reflected in Hpmulti and Hpinvmulti). See Forecast II, pp.69a-69b for details.
- if 1 % Liquidity effect prior on both MS and MD equations.
- [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior_covres_dobs(nvar,options_.ms.freq,options_.ms.nlags ,xdgel,mu,indxDummy,hpmsmd,indxmsmdeqn);
- else
- [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior(nvar,options_.ms.freq,options_.ms.nlags ,xdgel,mu);
- end
-
- %*** Combines asymmetric prior with linear restrictions
- [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Uiconst,Viconst,Pi,H0multi,Hpmulti,nvar);
-
- %*** Obtains the posterior matrices for estimation and inference
- [Pmat,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0invtld,Hpinvtld,Uiconst,Viconst);
-
- if options_.ms.contemp_reduced_form
- %*** Obtain the ML estimate
- A0hatinv = chol(H0inv{1}/fss); % upper triangular but lower triangular choleski
- A0hat=inv(A0hatinv);
- a0indx = find(A0hat);
- else
- %*** Obtain the ML estimate
- % load idenml
- x = 10*rand(sum(n0),1);
- H0 = eye(sum(n0));
- crit = 1.0e-9;
- nit = 10000;
- %
- [fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = ...
- csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Uiconst,nvar,n0,fss,H0inv);
-
- A0hat = fn_tran_b2a(xhat,Uiconst,nvar,n0)
- A0hatinv = inv(A0hat);
- fhat
- xhat
- grad
- itct
- fcount
- retcodehat
- save outm.mat xhat A0hat A0hatinv grad fhat itct itct fcount retcodehat
- end
- else
- %*** Obtain the posterior matrices for estimation and inference
- [Pmat,H0inv,Hpinv] = fn_dlrpostr(xtx,xty,yty,Uiconst,Viconst);
-
- if options_.ms.contemp_reduced_form
- %*** Obtain the ML estimate
- A0hatinv = chol(H0inv{1}/fss); % upper triangular but lower triangular choleski
- A0hat=inv(A0hatinv);
- a0indx = find(A0hat);
- else
- %*** Obtain the ML estimate
- % load idenml
- x = 10*rand(sum(n0),1);
- H0 = eye(sum(n0));
- crit = 1.0e-9;
- nit = 10000;
- %
- [fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = ...
- csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Uiconst,nvar,n0,fss,H0inv);
-
- A0hat = fn_tran_b2a(xhat,Uiconst,nvar,n0)
- A0hatinv = inv(A0hat);
- fhat
- xhat
- grad
- itct
- fcount
- retcodehat
- save outm.mat xhat A0hat A0hatinv grad fhat itct itct fcount retcodehat
- end
- end
-
- %**** impulse responses
- swish = A0hatinv; % each column corresponds to an equation
- if options_.ms.contemp_reduced_form
- xhat = A0hat(a0indx);
- Bhat=Pmat{1};
- Fhat = Bhat*A0hat
- ghat = NaN;
- else
- xhat = fn_tran_a2b(A0hat,Uiconst,nvar,n0);
- [Fhat,ghat] = fn_gfmean(xhat,Pmat,Viconst,nvar,ncoef,n0,np);
- if options_.ms.cross_restrictions
- Fhatur0P = Fhat; % ur: unrestriced across A0 and A+
- for ki = 1:size(ixmC0Pres,1) % loop through the number of equations in which
- % cross-A0-A+ restrictions occur. See St. Louis Note p.5.
- ixeq = ixmC0Pres{ki}(1,1); % index for the jth equation in consideration.
- Lit = Viconst{ixeq}(ixmC0Pres{ki}(:,2),:); % transposed restriction matrix Li
- % V_j(i,:) in f_j(i) = V_j(i,:)*g_j
- ci = ixmC0Pres{ki}(:,4) .* A0hat(ixmC0Pres{ki}(:,3),ixeq);
- % s * a_j(h) in the restriction f_j(i) = s * a_j(h).
- LtH = Lit/Hpinv{ixeq};
- HLV = LtH'/(LtH*Lit');
- gihat = Viconst{ixeq}'*Fhatur0P(:,ixeq);
- Fhat(:,ixeq) = Viconst{ixeq}*(gihat + HLV*(ci-Lit*gihat));
- end
- end
- Fhat
- Bhat = Fhat/A0hat; % ncoef-by-nvar reduced form lagged parameters.
- end
- nn = [nvar options_.ms.nlags imstp];
- imfhat = fn_impulse(Bhat,swish,nn); % in the form that is congenial to RATS
- imf3hat=reshape(imfhat,size(imfhat,1),nvar,nvar);
- % imf3: row--steps, column--nvar responses, 3rd dimension--nvar shocks
- imf3shat=permute(imf3hat,[1 3 2]);
- % imf3s: permuted so that row--steps, column--nvar shocks,
- % 3rd dimension--nvar responses
- % Note: reshape(imf3s(1,:,:),nvar,nvar) = A0in (columns -- equations)
- if options_.ms.indxgimfhat
- figure
- end
- scaleout = fn_imcgraph(imfhat,nvar,imstp,xlab,ylab,options_.ms.indxgimfhat);
- imfstd = max(abs(scaleout)'); % row: nvar (largest number); used for standard deviations
-
- %
- % %**** save stds. of both data and impulse responses in idfile1
- % temp = [std(yactqmyge(:,3:end)); std(yactyrge(:,3:end)); imfstd]; %<<>>
- % save idenyimstd.prn temp -ascii % export forecast and impulse response to the file "idenyimstd.prn", 3-by-nvar
- % %
- % %**** save stds. of both data and impulse responses in idfile1
- % temp = [std(yactqmyge(:,3:end)); std(yactyrge(:,3:end)); imfstd]; %<<>>
- % save idenyimstd.prn temp -ascii % export forecast and impulse response to the file "idenyimstd.prn", 3-by-nvar
- % if options_.ms.indxparr
- % idfile1='idenyimstd';
- % end
-
- %=====================================
- % Now, out-of-sample forecasts. Note: Hm1t does not change with A0.
- %=====================================
- %
- % * updating the last row of X (phi) with the current (last row of) y.
- tcwx = nvar*options_.ms.nlags ; % total coefficeint without exogenous variables
- phil = phi(size(phi,1),:);
- phil(nvar+1:tcwx) = phil(1:tcwx-nvar);
- phil(1:nvar) = y(end,:);
- %*** exogenous variables excluding constant terms
- if (nexo>1)
- Xexoe = fn_dataext([yrEnd qmEnd],[yrEnd qmEnd],xdatae(:,[1:2 2+nvar+1:2+nvar+nexo-1]));
- phil(1,tcwx+1:tcwx+nexo-1) = Xexoe(1,3:end);
- end
- %
- %*** ML unconditional point forecast
- nn = [nvar options_.ms.nlags nfqm];
- if nexo<2
- yforehat = fn_forecast(Bhat,phil,nn); % nfqm-by-nvar, in log
- else
- Xfexoe = fn_dataext(fdates(1,:),fdates(numel(fdates),:),xdatae(:,[1:2 2+nvar+1:2+nvar+nexo-1]));
- %Xfexoe = fn_dataext(fdates(1,:),fdates(end,:),xdatae(:,[1:2 2+nvar+1:2+nvar+nexo-1]));
- yforehat = fn_forecast(Bhat,phil,nn,nexo,Xfexoe(:,3:end)); % nfqm-by-nvar, in log
- end
- yforehate = [fdates yforehat];
- %
- yact1e = fn_dataext([yrEnd-nayr 1],[yrEnd qmEnd],xdatae(:,1:nvar+2));
- if options_.ms.real_pseudo_forecast
- %yact2e = fn_dataext([yrEnd-nayr 1],E2yrqm,xdatae);
- yact2e = fn_dataext([yrEnd-nayr 1],[fdates(end,1) options_.ms.freq],xdatae(:,1:nvar+2));
- else
- yact2e=yact1e;
- end
- yafhate = [yact1e; yforehate]; % actual and forecast
- %
- %===== Converted to mg, qg, and calendar yg
- %
- [yafyrghate,yafyrhate,yafqmyghate] = fn_datana(yafhate,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno));
- % actual and forecast growth rates
- [yact2yrge,yact2yre,yact2qmyge] = fn_datana(yact2e,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno));
- % only actual growth rates
- yafyrghate
- if options_.ms.indxgforhat
- keyindx = [1:nvar];
- conlab=['unconditional'];
-
- figure
- yafyrghate(:,3:end) = yafyrghate(:,3:end)/100;
- yact2yrge(:,3:end) = yact2yrge(:,3:end)/100;
- fn_foregraph(yafyrghate,yact2yrge,keyindx,rnum,cnum,options_.ms.freq,ylab,forelabel,conlab)
- end
-
- %-------------------------------------------------
- % Setup for point conditional forecast
- % ML Conditional Forecast
- %-------------------------------------------------
- %
- %% See Zha's note "Forecast (1)" p. 5, RATS manual (some errors in RATS), etc.
- %
- %% Some notations: y(t+1) = y(t)B1 + e(t+1)inv(A0). e(t+1) is 1-by-n.
- %% Let r(t+1)=e(t+1)inv(A0) + e(t+2)C + .... where inv(A0) is impulse
- %% response at t=1, C at t=2, etc. The row of inv(A0) or C is
- %% all responses to one shock.
- %% Let r be q-by-1 (such as r(1) = r(t+1)
- %% = y(t+1) (constrained) - y(t+1) (forecast)).
- %% Use impulse responses to find out R (k-by-q) where k=nvar*nsteps
- %% where nsteps the largest constrained step. The key of the program
- %% is to creat R using impulse responses
- %% Optimal solution for shock e where R'*e=r and e is k-by-1 is
- %% e = R*inv(R'*R)*r.
- %
-
- if (nconstr > 0)
- %*** initializing
- stepcon=cell(nconstr,1); % initializing, value y conditioned
- valuecon=zeros(nconstr,1); % initializing, value y conditioned
- varcon=zeros(nconstr,1); % initializing, endogous variables conditioned
- varcon(:)=PorR; % 1: Pcm; 3: FFR; 5: CPI
-
- %
- for i=1:nconstr
- if i<=nconstr1
- stepcon{i}=i; % FFR
- elseif i<=2*nconstr1
- stepcon{i}=i-nconstr1; % FFR
- elseif i<=3*nconstr1
- stepcon{i}=i-2*nconstr1; % FFR
- elseif i<=4*nconstr1
- stepcon{i}=i-3*nconstr1; % FFR
- elseif i<=5*nconstr1
- stepcon{i}=i-4*nconstr1; % FFR
- elseif i<=6*nconstr1
- stepcon{i}=i-5*nconstr1; % FFR
- end
- end
-
-% for i=1:nconstr
-% stepcon{i}=i; % FFR
-% end
-
-% bend=12;
-% stepcon{1}=[1:bend]'; % average over
-% stepcon{nconstr1+1}=[1:options_.ms.freq-qmSub]'; % average over the remaing months in 1st forecast year
-% stepcon{nconstr1+2}=[options_.ms.freq-qmSub+1:options_.ms.freq-qmSub+12]'; % average over 12 months next year
-% stepcon{nconstr1+3}=[options_.ms.freq-qmSub+13:options_.ms.freq-qmSub+24]'; % average over 12 months. 3rd year
-% stepcon{nconstr1+4}=[options_.ms.freq-qmSub+25:options_.ms.freq-qmSub+36]'; % average over 12 months. 4th year
-
-% %**** avearage condition over, say, options_.ms.freq periods
-% if qmEnd==options_.ms.freq
-% stepcon{1}=[1:options_.ms.freq]'; % average over the remaing periods in 1st forecast year
-% else
-% stepcon{1}=[1:options_.ms.freq-qmEnd]'; % average over the remaing periods in 1st forecast year
-% end
-% for kj=2:nconstr
-% stepcon{kj}=[length(stepcon{kj-1})+1:length(stepcon{kj-1})+options_.ms.freq]'; % average over 12 months next year
-% end
-
- if options_.ms.real_pseudo_forecast
-% %*** conditions in every period
-% for i=1:nconstr
-% valuecon(i) = yact(actup+i,varcon(i));
-% %valuecon(i) = mean( yact(actup+1:actup+bend,varcon(i)) );
-% %valuecon(i) = 0.060; % 95:01
-% %valuecon(i) = (0.0475+0.055)/2; % 94:10
-% end
-
-% %*** average condtions over,say, options_.ms.freq periods.
-% for i=nconstr1+1:nconstr1+nconstr2
-% i=1;
-% valuecon(nconstr1+i) = ( ( mean(ylast12Cal(:,varcon(nconstr1+i)),1) + ...
-% log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100) )*options_.ms.freq - ...
-% yCal_1(:,varcon(nconstr1+i)) ) ./ length(stepcon{nconstr1+i});
-% % the same as unconditional "yactCalyg" 1st calendar year
-% i=2;
-% valuecon(nconstr1+i) = mean(ylast12Cal(:,varcon(nconstr1+i))) + ...
-% log(1+yactCalyg(yAg-yFg+1,varcon(nconstr1+i))/100) ...
-% + log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100);
-% % the same as actual "yactCalgy" 2nd calendar year
-% i=3;
-% valuecon(nconstr1+i) = valuecon(nconstr1+i-1) + ...
-% log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100);
-% % the same as actual "yactCalgy" 3rd calendar year
-% %i=4;
-% %valuecon(nconstr1+i) = valuecon(nconstr1+i-1) + ...
-% % log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100);
-% % the same as actual "yactCalgy" 4th calendar year
-% end
-
- %*** conditions in every period
- vpntM = fn_dataext(E1yrqm, E2yrqm,xdatae); % point value matrix with dates
- % vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates
- for i=1:nconstr
- if i<=nconstr1
- valuecon(i) = vpntM(i,2+varcon(i)); % 2: first 2 elements are dates
- elseif i<=2*nconstr1
- valuecon(i) = vpntM(i-nconstr1,2+varcon(i));
- elseif i<=3*nconstr1
- valuecon(i) = vpntM(i-2*nconstr1,2+varcon(i));
- elseif i<=4*nconstr1
- valuecon(i) = vpntM(i-3*nconstr1,2+varcon(i));
- elseif i<=5*nconstr1
- valuecon(i) = vpntM(i-4*nconstr1,2+varcon(i));
- elseif i<=6*nconstr1
- valuecon(i) = vpntM(i-5*nconstr1,2+varcon(i));
- end
- end
-
-% %*** average condtions over,say, options_.ms.freq periods.
-% if qmEnd==options_.ms.freq
-% vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates
-% valuecon(1) = vaveM(1,2+varcon(1)); % 2: first 2 elements are dates
-% else
-% vaveM = fn_dataext([yrEnd 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates
-% yactrem = fn_dataext([yrEnd qmEnd+1],[yrEnd options_.ms.freq],xdatae);
-% valuecon(1) = sum(yactrem(:,2+varcon(1)),1)/length(stepcon{1});
-% % 2: first 2 elements are dates
-% end
-% for kj=2:nconstr
-% valuecon(kj) = vaveM(kj,2+varcon(kj)); % 2: first 2 elements are dates
-% end
- else
- vpntM = dataext([yrEnd qmEnd+1],[yrEnd qmEnd+2],xdatae); % point value matrix with dates
- for i=1:nconstr
- if i<=nconstr1
- valuecon(i) = vpntM(i,2+varcon(i)); % 2: first 2 elements are dates; Poil
- elseif i<=2*nconstr1
- valuecon(i) = vpntM(i-nconstr1,2+varcon(i)); % 2: first 2 elements are dates; M2
- elseif i<=3*nconstr1
- valuecon(i) = vpntM(i-2*nconstr1,2+varcon(i)); % 2: first 2 elements are dates; FFR
- elseif i<=4*nconstr1
- valuecon(i) = vpntM(i-3*nconstr1,2+varcon(i)); % 2: first 2 elements are dates; CPI
- elseif i<=5*nconstr1
- valuecon(i) = vpntM(i-4*nconstr1,2+varcon(i)); % 2: first 2 elements are dates; U
- elseif i<=5*nconstr1+nconstr2
- valuecon(i)=xdata(end,5)+(i-5*nconstr1)*log(1.001)/options_.ms.freq; %CPI
- elseif i<=5*nconstr1+2*nconstr2
- valuecon(i)=0.0725; %FFR
- else
- valuecon(i)=xdata(end,6)+(i-5*nconstr1-2*nconstr2)*0.01/nfqm; %U
- end
- end
- %valuecon(i) = 0.060; % 95:01
- end
- else
- valuecon = [];
- stepcon = [];
- varcon = [];
- end
-
- nstepsm = 0; % initializing, the maximum step in all constraints
- for i=1:nconstr
- nstepsm = max([nstepsm max(stepcon{i})]);
- end
-
- if cnum
- if options_.ms.real_pseudo_forecast && options_.ms.banact
- for i=1:length(banstp{1})
- banval{1}(1:length(banstp{1}),1) = ...
- yactCalyg(yAg-yFg+1:yAg-yFg+length(banstp{1}),banvar(1)) - 2;
- banval{1}(1:length(banstp{1}),2) = ...
- yactCalyg(yAg-yFg+1:yAg-yFg+length(banstp{1}),banvar(1)) + 2;
- end
- end
- end
-
-
- %===================================================
- % ML conditional forecast
- %===================================================
- %/*
- [ychat,Estr,rcon] = fn_fcstidcnd(valuecon,stepcon,varcon,nstepsm,...
- nconstr,options_.ms.eq_ms,nvar,options_.ms.nlags ,phil,0,0,yforehat,imf3shat,A0hat,Bhat,...
- nfqm,options_.ms.tlindx,options_.ms.tlnumber,options_.ms.ncms,options_.ms.eq_cms);
- ychate = [fdates ychat];
- yachate = [yact1e; ychate]; % actual and condtional forecast
- %===== Converted to mg, qg, and calendar yg
- [yacyrghate,yacyrhate,yacqmyghate] = fn_datana(yachate,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno));
- % actual and conditional forecast growth rates
- if options_.ms.indxgdls && nconstr
- keyindx = [1:nvar];
- % conlab=['conditional on' ylab{PorR(1)}];
- conlab=['v-conditions'];
-
- figure
- fn_foregraph(yafyrghate,yact2yrge,keyindx,rnum,cnum,options_.ms.freq,ylab,forelabel,conlab)
- end
-
- if options_.ms.ncsk
- Estr = zeros(nfqm,nvar);
- Estr(1:2,:) = [
- -2.1838 -1.5779 0.53064 -0.099425 -0.69269 -1.0391
- 1.9407 3.3138 -0.10563 -0.55457 -0.68772 1.3534
- ];
- Estr(3:6,3) = [0.5*ones(1,4)]'; % MD shocks
-
- Estr(3:10,2) = [1.5 1.5 1.5*ones(1,6)]'; % MS shocks
-
- %Estr(3:6,6) = 1*ones(4,1); % U shocks
- %Estr(8:11,4) = 1*ones(4,1); % y shocks
-
- %Estr(3:10,2) = [2.5 2.5 1.5*ones(1,6)]'; % MS shocks alone
-
- nn = [nvar options_.ms.noptions_.ms.nlags nfqm];
- ycEhat = forefixe(A0hat,Bhat,phil,nn,Estr);
- ycEhate = [fdates ycEhat];
- yacEhate = [yact1e; ycEhate]; % actual and condtional forecast
- %===== Converted to mg, qg, and calendar yg
- [yacEyrghate,yacEyrhate,yacEqmyghate] = datana(yacEhate,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno));
- % actual and conditional forecast growth rates
- disp([sprintf('%4.0f %2.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',yacEyrghate')])
-
- if 1
- keyindx = [1:nvar];
- % conlab=['conditional on' ylab{PorR(1)}];
- conlab=['shock-conditions'];
-
- figure
- gyrfore(yacEyrghate,yact2yrge,keyindx,rnum,cnum,ylab,forelabel,conlab)
- end
- end
-
- %-----------------------------------------------------------
- % Compute structural shocks for the whole sample period excluding dummy observations.
- %-----------------------------------------------------------
- ywod = y(options_.ms.dummy_obs +1:end,:); % without dummy observations
- phiwod=phi(options_.ms.dummy_obs +1:end,:); % without dummy observations
- eplhat=ywod*A0hat-phiwod*Fhat;
- qmStartWod = mod(qmStart+options_.ms.nlags ,options_.ms.freq);
- if (~qmStartWod)
- qmStartWod = options_.ms.freq;
- end
- yrStartWod = yrStart + floor((qmStart+options_.ms.nlags -1)/options_.ms.freq);
- dateswod = fn_dataext([yrStartWod qmStartWod],[yrEnd qmEnd],xdatae(:,[1:2]));
- eplhate = [dateswod eplhat];
-
- Aphat = Fhat;
-end
-
-%----------------------------------------
-% Tests for LR, HQ, Akaike, Schwarz. The following gives a guidance.
-% But the computation has to be done in a different M file by exporting fhat's
-% from different idfile's.
-%----------------------------------------
-%
-%if ~options_.ms.contemp_reduced_form
-% SpHR=A0in'*A0in;
-%end
-%%
-%if ~isnan(SpHR) && ~options_.ms.contemp_reduced_form
-% warning(' ')
-% disp('Make sure you run the program with options_.ms.contemp_reduced_form =1 first.')
-% disp('Otherwise, the following test results such as Schwartz are incorrect.')
-% disp('All other results such as A0ml and imfs, however, are correct.')
-% disp('Press anykey to contintue or ctrl-c to stop now')
-% pause
-
-% load SpHUout
-
-% logLHU=-fss*sum(log(diag(chol(SpHU)))) -0.5*fss*nvar % unrestricted logLH
-
-% logLHR=-fhat % restricted logLH
-% tra = reshape(SpHU,nvar*nvar,1)'*reshape(A0*A0',nvar*nvar,1);
-% df=(nvar*(nvar+1)/2 - length(a0indx));
-% S=2*(logLHU-logLHR);
-% SC = (nvar*(nvar+1)/2 - length(a0indx)) * log(fss);
-% disp(['T -- effective sample size: ' num2str(fss)])
-% disp(['Trace in the overidentified posterior: ' num2str(tra)])
-% disp(['Chi2 term -- 2*(logLHU-logLHR): ' num2str(S)])
-% disp(['Degrees of freedom: ' num2str(df)])
-% disp(['SC -- df*log(T): ' num2str(SC)])
-% disp(['Akaike -- 2*df: ' num2str(2*df)])
-% disp(['Classical Asymptotic Prob at chi2 term: ' num2str(cdf('chi2',S,df))])
-
-% %*** The following is the eigenanalysis in the difference between
-% %*** unrestricted (U) and restricted (R)
-% norm(A0'*SpHU*A0-diag(diag(ones(6))))/6;
-% norm(SpHU-A0in'*A0in)/6;
-
-% corU = corr(SpHU);
-% corR = corr(SpHR);
-
-% [vU,dU]=eigsort(SpHU,1);
-% [vR,dR]=eigsort(SpHR,1);
-
-% [log(diag(dU)) log(diag(dR)) log(diag(dU))-log(diag(dR))];
-
-% sum(log(diag(dU)));
-% sum(log(diag(dR)));
-%else
-% disp('To run SC test, turn options_.ms.contemp_reduced_form =1 first and then turn options_.ms.contemp_reduced_form =0')
-%end
-
-
-%***** Simply regression
-%X=[phi(:,3) y(:,2)-phi(:,2) y(:,1)-phi(:,7) ones(fss,1)];
-%� Y=y(:,3);
-%� b=regress(Y,X)
-
-%=== Computes the roots for the whole system.
-rootsinv = fn_varoots(Bhat,nvar,options_.ms.nlags )
-abs(rootsinv)
-
-
-bhat =xhat;
-n0const=n0; % For constant parameter models.
-n0const=n0; % For constant parameter models.
-npconst=np; % For constant parameter models.
-save outdata_a0dp_const.mat A0hat bhat Aphat n0const
+%function msstart2(options_)
+% This .m file is called for point graphics or error bands and
+% It starts for both data and Bayesian estimation (when IxEstimat==0,
+% no estimation but only data analysis), which allows you to set
+% (a) available data range,
+% (b) sample range,
+% (c) rearrangement of actual data such as mlog, qg, yg
+% (d) conditions of shocks 'Cms',
+% (c) conditions of variables 'nconstr',
+% (e) soft conditions 'nbancon,'
+% (f) produce point conditional forecast (at least conditions on variables).
+%
+% February 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 .
+
+% ** ONLY UNDER UNIX SYSTEM
+%path(path,'/usr2/f1taz14/mymatlab')
+
+%addpath('C:\SoftWDisk\MATLAB6p5\toolbox\cstz')
+
+
+msstart_setup
+%===========================================
+% Exordium II
+%===========================================
+%options_.ms.ncsk = 0; % conditional directly on shoocks. Unlike Cms, not on variables that back out shocks
+%options_.ms.nstd = 6; % number of standard deviations to cover the range in which distributions are put to bins
+%options_.ms.ninv = 1000; % the number of bins for grouping, say, impulse responses
+%options_.ms.Indxcol = [1:nvar]; % a vector of random columns in which MC draws are made.
+%
+%options_.ms.indxparr = 1; % 1: parameters random; 0: no randomness in parameters
+ % Note, when 0, there is no effect from the values of options_.ms.IndxAp, options_.ms.Aband, etc.
+%options_.ms.indxovr = 0; % 1: distributions for other variables of interest; 0: no distribution.
+ % Example: joint distribution of a(1) and a(2). Only for specific purposes
+%options_.ms.Aband = 1; % 1: error bands with only A0 and A+ random.
+%options_.ms.IndxAp = 1; % 1: generate draws of A+; 0: no such draws.
+ % Note: when options_.ms.IndxAp=0, there is no effect from the values of options_.ms.options_.ms.options_.ms.options_.ms.indximf, IndxFore,
+ % or options_.ms.apband.
+%options_.ms.apband = 1; % 1: error bands for A+; 0: no error bands for A+.
+%*** The following (impulse responses and forecasts) is used only if options_.ms.IndxAp=1
+%options_.ms.indximf = 1; % 1: generate draws of impulse responses; 0: no such draws (thus no effect
+ % from options_.ms.imfband)
+%options_.ms.imfband = 1; % 1: error bands for impulse responses; 0: no error bands
+%options_.ms.indxfore = 0; % 1: generate draws of forecasts; 0: no such draws (thus no effect from options_.ms.foreband)
+%options_.ms.foreband = 0; % 1: error bands for out-of-sample forecasts; 0: no error bands
+%
+%options_.ms.indxgforhat = 1; % 1: plot unconditoinal forecasts; 0: no such plot
+rnum = nvar; % number of rows in the graph
+cnum = 1; % number of columns in the graph
+if rnum*cnum>
+nconstr=0 ; %6*nconstr1;
+options_.ms.eq_ms = []; % location of MS equation; if [], all shocks
+PorR = [4*ones(nconstr1,1);2*ones(nconstr1,1);3*ones(nconstr1,1)]; % the variable conditioned. 1: Pcm; 3: FFR; 4: CPI
+PorR = [PorR;1*ones(nconstr1,1);5*ones(nconstr1,1);6*ones(nconstr1,1)];
+%PorR = [3 5]; % the variable conditioned. 3: FFR; 5: CPI
+%PorR = 3;
+
+%%%%----------------------------------------
+% Conditions directly on future shocks
+%
+%options_.ms.cms = 0 % 1: condition on ms shocks; 0: disable this and "fidcnderr.m" gives
+ % unconditional forecasts if nconstr = 0 as well; <<>>
+%options_.ms.ncms = 0; % number of the stance of policy; 0 if no tightening or loosening
+%options_.ms.eq_cms = 1; % location of MS shocks
+options_.ms.tlindx = 1*ones(1,options_.ms.ncms); % 1-by-options_.ms.ncms vector; 1: tightening; 0: loosen
+options_.ms.tlnumber = [0.5 0.5 0 0]; %94:4 % [2 2 1.5 1.5]; %79:9 %[1.5 1.5 1 1]; 90:9
+ % 1-by-options_.ms.ncms vector; cut-off point for MS shocks
+TLmean = zeros(1,options_.ms.ncms);
+ % unconditional, i.e., 0 mean, for the final report in the paper
+if options_.ms.cms
+ options_.ms.eq_ms = [];
+ % At least at this point, it makes no sense to have DLS type of options_.ms.eq_ms; 10/12/98
+ if all(isfinite(options_.ms.tlnumber))
+ for k=1:options_.ms.ncms
+ TLmean(k) = lcnmean(options_.ms.tlnumber(k),options_.ms.tlindx(k));
+ % shock mean magnitude. 1: tight; 0: loose
+ % Never used for any subsequent computation but
+ % simply used for the final report in the paper.
+ %options_.ms.tlnumber(k) = fzero('lcutoff',0,[],[],TLmean(k))
+ % get an idea about the cutoff point given TLmean instead
+
+ end
+ end
+else
+ options_.ms.ncms = 0; % only for the use of the graph by msprobg.m
+ options_.ms.tlnumber = NaN*ones(1,options_.ms.ncms);
+ % -infinity, only for the use of the graph by msprobg.m
+end
+
+
+%%%%----------------------------------------
+% Soft conditions on variables
+%
+%cnum = 0 % # of band condtions; when 0, disable this option
+ % Note (different from "fidencon") that each condition corres. to variable
+%options_.ms.banact = 1; % 1: use infor on actual; 0: preset without infor on actual
+if cnum
+ banindx = cell(cnum,1); % index for each variable or conditon
+ banstp = cell(cnum,1); % steps: annual in general
+ banvar = zeros(cnum,1); % varables: annual in general
+ banval = cell(cnum,1); % band value (each variable occupy a cell)
+ badval{1} = zeros(length(banstp{1}),2); % 2: lower or higher bound
+
+ banstp{1} = 1:4; % 3 or 4 years
+ banvar(1) = 3; % 3: FFR; 5: CPI
+ if ~options_.ms.banact
+ for i=1:length(banstp{1})
+ banval{1}(i,:) = [5.0 10.0];
+ end
+ end
+end
+%
+pause(1)
+disp(' ')
+disp('For uncondtional forecasts, set nconstr = options_.ms.cms = cnum = 0')
+pause(1)
+%
+%=================================================
+%====== End of exordium ==========================
+%=================================================
+
+
+
+
+
+%(1)--------------------------------------
+% Further data analysis
+%(1)--------------------------------------
+%
+if (options_.ms.freq==12)
+ nStart=(yrStart-options_.ms.initial_year )*12+qmStart-options_.ms.initial_subperiod ; % positive number of months at the start
+ nEnd=(yrEnd-options_.ms.final_year )*12+qmEnd-options_.ms.final_subperiod ; % negative number of months towards end
+elseif (options_.ms.freq==4)
+ nStart=(yrStart-options_.ms.initial_year )*4+qmStart-options_.ms.initial_subperiod ; % positive number of months at the start
+ nEnd=(yrEnd-options_.ms.final_year )*4+qmEnd-options_.ms.final_subperiod ; % negative number of months towards end
+else
+ error('Error: this code is only good for monthly/quarterly data!!!')
+ return
+end
+%
+if nEnd>0 || nStart<0
+ disp('Warning: this particular sample consider is out of bounds of the data!!!')
+ return
+end
+%*** Note, both xdgel and xdata have the same start with the specific sample
+xdgel=options_.data(nStart+1:nData+nEnd,options_.ms.vlist);
+ % gel: general options_.data within sample (nSample)
+if ~(nSample==size(xdgel,1))
+ warning('The sample size (including options_.ms.nlags ) and data are incompatible')
+ disp('Check to make sure nSample and size(xdgel,1) are the same')
+ return
+end
+%
+baddata = find(isnan(xdgel));
+if ~isempty(baddata)
+ warning('Some data for this selected sample are actually unavailable.')
+ disp('Hit any key to continue, or ctrl-c to abort')
+ pause
+end
+%
+if options_.ms.initial_subperiod ==1
+ yrB = options_.ms.initial_year ; qmB = options_.ms.initial_subperiod ;
+else
+ yrB = options_.ms.initial_year +1; qmB = 1;
+end
+yrF = options_.ms.final_year ; qmF = options_.ms.final_subperiod ;
+[Mdate,tmp] = fn_calyrqm(options_.ms.freq,[options_.ms.initial_year options_.ms.initial_subperiod ],[options_.ms.final_year options_.ms.final_subperiod ]);
+xdatae=[Mdate options_.data(1:nData,options_.ms.vlist)];
+ % beyond sample into forecast horizon until the end of the data options_.ms.final_year :options_.ms.final_subperiod
+ % Note: may contain NaN data. So must be careful about its use
+
+%=========== Obtain prior-period, period-to-last period, and annual growth rates
+[yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(xdatae,options_.ms.freq,options_.ms.log_var,options_.ms.percent_var,[yrB qmB],[yrF qmF]);
+qdates = zeros(size(yactqmyge,1),1);
+for ki=1:length(qdates)
+ qdates(ki) = yactqmyge(1,1) + (yactqmyge(1,2)+ki-2)/options_.ms.freq;
+end
+for ki=1:nvar
+ figure
+ plot(qdates, yactqmyge(:,2+ki)/100)
+ xlabel(options_.ms.varlist{ki})
+end
+save outactqmygdata.prn yactqmyge -ascii
+
+
+
+%=========== Write the output on the screen or to a file in an organized way ==============
+%disp([sprintf('%4.0f %2.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',yactyrge')])
+spstr1 = 'disp([sprintf(';
+spstr2 = '%4.0f %2.0f';
+yactyrget=yactyrge';
+for ki=1:length(options_.ms.vlist)
+ if ki==length(options_.ms.vlist)
+ spstr2 = [spstr2 ' %8.3f\n'];
+ else
+ spstr2 = [spstr2 ' %8.3f'];
+ end
+end
+spstr = [spstr1 'spstr2' ', yactyrget)])'];
+eval(spstr)
+
+%
+fid = fopen('outyrqm.prn','w');
+%fprintf(fid,'%4.0f %2.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',yactyrge');
+fpstr1 = 'fprintf(fid,';
+fpstr2 = '%4.0f %2.0f';
+for ki=1:nvar
+ if ki==nvar
+ fpstr2 = [fpstr2 ' %8.3f\n'];
+ else
+ fpstr2 = [fpstr2 ' %8.3f'];
+ end
+end
+fpstr = [fpstr1 'fpstr2' ', yactyrget);'];
+eval(fpstr)
+fclose(fid);
+
+
+
+if options_.ms.indxestima
+ %(2)----------------------------------------------------------------------------
+ % Estimation
+ % ML forecast and impulse responses
+ % Hard or soft conditions for conditional forecasts
+ %(2)----------------------------------------------------------------------------
+ %
+ %* Arranged data information, WITHOUT dummy obs when 0 after mu is used. See fn_rnrprior_covres_dobs.m for using the dummy
+ % observations as part of an explicit prior.
+ [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh] = fn_dataxy(nvar,options_.ms.nlags ,xdgel,mu,0,nexo);
+ if qmStart+options_.ms.nlags -options_.ms.dummy_obs >0
+ qmStartEsti = rem(qmStart+options_.ms.nlags -options_.ms.dummy_obs ,options_.ms.freq); % dummy observations are included in the sample.
+ if (~qmStartEsti)
+ qmStartEsti = options_.ms.freq;
+ end
+ yrStartEsti = yrStart + floor((qmStart+options_.ms.nlags -options_.ms.dummy_obs )/(options_.ms.freq+0.01));
+ % + 0.01 (or any number < 1) is used so that qmStart+options_.ms.nlags -options_.ms.dummy_obs ==?*options_.ms.freq doesn't give us an extra year forward.
+ else
+ qmStartEsti = options_.ms.freq + rem(qmStart+options_.ms.nlags -options_.ms.dummy_obs ,options_.ms.freq); % dummy observations are included in the sample.
+ if (qmStart+options_.ms.nlags -options_.ms.dummy_obs ==0)
+ yrStartEsti = yrStart - 1; % one year back.
+ else
+ yrStartEsti = yrStart + floor((qmStart+options_.ms.nlags -options_.ms.dummy_obs )/(options_.ms.freq-0.01));
+ % - 0.01 (or any number < 1) is used so that qmStart+options_.ms.nlags -options_.ms.dummy_obs ==-?*options_.ms.freq give us an extra year back.
+ end
+ end
+ dateswd = fn_dataext([yrStartEsti qmStartEsti],[yrEnd qmEnd],xdatae(:,[1:2])); % dates with dummies
+ phie = [dateswd phi];
+ ye = [dateswd y];
+
+ %* Obtain linear restrictions
+ [Uiconst,Viconst,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,nvar,nexo,options_.ms );
+ if min(n0)==0
+ disp(' ')
+ warning('A0: restrictions in dlrprior.m give no free parameter in one of equations')
+ disp('Press ctrl-c to abort')
+ pause
+ elseif min(np)==0
+ disp(' ')
+ warning('Ap: Restrictions in dlrprior.m give no free parameter in one of equations')
+ disp('Press ctrl-c to abort')
+ pause
+ end
+
+ if options_.ms.contemp_reduced_form
+ Uiconst=cell(nvar,1); Viconst=cell(ncoef,1);
+ for kj=1:nvar
+ Uiconst{kj} = eye(nvar); Viconst{kj} = eye(ncoef);
+ end
+ end
+
+ if options_.ms.bayesian_prior
+ %*** Obtains asymmetric prior (with no linear restrictions) with dummy observations as part of an explicit prior (i.e,
+ % reflected in Hpmulti and Hpinvmulti). See Forecast II, pp.69a-69b for details.
+ if 1 % Liquidity effect prior on both MS and MD equations.
+ [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior_covres_dobs(nvar,options_.ms.freq,options_.ms.nlags ,xdgel,mu,indxDummy,hpmsmd,indxmsmdeqn);
+ else
+ [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior(nvar,options_.ms.freq,options_.ms.nlags ,xdgel,mu);
+ end
+
+ %*** Combines asymmetric prior with linear restrictions
+ [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Uiconst,Viconst,Pi,H0multi,Hpmulti,nvar);
+
+ %*** Obtains the posterior matrices for estimation and inference
+ [Pmat,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0invtld,Hpinvtld,Uiconst,Viconst);
+
+ if options_.ms.contemp_reduced_form
+ %*** Obtain the ML estimate
+ A0hatinv = chol(H0inv{1}/fss); % upper triangular but lower triangular choleski
+ A0hat=inv(A0hatinv);
+ a0indx = find(A0hat);
+ else
+ %*** Obtain the ML estimate
+ % load idenml
+ x = 10*rand(sum(n0),1);
+ H0 = eye(sum(n0));
+ crit = 1.0e-9;
+ nit = 10000;
+ %
+ [fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = ...
+ csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Uiconst,nvar,n0,fss,H0inv);
+
+ A0hat = fn_tran_b2a(xhat,Uiconst,nvar,n0)
+ A0hatinv = inv(A0hat);
+ fhat
+ xhat
+ grad
+ itct
+ fcount
+ retcodehat
+ save outm.mat xhat A0hat A0hatinv grad fhat itct itct fcount retcodehat
+ end
+ else
+ %*** Obtain the posterior matrices for estimation and inference
+ [Pmat,H0inv,Hpinv] = fn_dlrpostr(xtx,xty,yty,Uiconst,Viconst);
+
+ if options_.ms.contemp_reduced_form
+ %*** Obtain the ML estimate
+ A0hatinv = chol(H0inv{1}/fss); % upper triangular but lower triangular choleski
+ A0hat=inv(A0hatinv);
+ a0indx = find(A0hat);
+ else
+ %*** Obtain the ML estimate
+ % load idenml
+ x = 10*rand(sum(n0),1);
+ H0 = eye(sum(n0));
+ crit = 1.0e-9;
+ nit = 10000;
+ %
+ [fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = ...
+ csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Uiconst,nvar,n0,fss,H0inv);
+
+ A0hat = fn_tran_b2a(xhat,Uiconst,nvar,n0)
+ A0hatinv = inv(A0hat);
+ fhat
+ xhat
+ grad
+ itct
+ fcount
+ retcodehat
+ save outm.mat xhat A0hat A0hatinv grad fhat itct itct fcount retcodehat
+ end
+ end
+
+ %**** impulse responses
+ swish = A0hatinv; % each column corresponds to an equation
+ if options_.ms.contemp_reduced_form
+ xhat = A0hat(a0indx);
+ Bhat=Pmat{1};
+ Fhat = Bhat*A0hat
+ ghat = NaN;
+ else
+ xhat = fn_tran_a2b(A0hat,Uiconst,nvar,n0);
+ [Fhat,ghat] = fn_gfmean(xhat,Pmat,Viconst,nvar,ncoef,n0,np);
+ if options_.ms.cross_restrictions
+ Fhatur0P = Fhat; % ur: unrestriced across A0 and A+
+ for ki = 1:size(ixmC0Pres,1) % loop through the number of equations in which
+ % cross-A0-A+ restrictions occur. See St. Louis Note p.5.
+ ixeq = ixmC0Pres{ki}(1,1); % index for the jth equation in consideration.
+ Lit = Viconst{ixeq}(ixmC0Pres{ki}(:,2),:); % transposed restriction matrix Li
+ % V_j(i,:) in f_j(i) = V_j(i,:)*g_j
+ ci = ixmC0Pres{ki}(:,4) .* A0hat(ixmC0Pres{ki}(:,3),ixeq);
+ % s * a_j(h) in the restriction f_j(i) = s * a_j(h).
+ LtH = Lit/Hpinv{ixeq};
+ HLV = LtH'/(LtH*Lit');
+ gihat = Viconst{ixeq}'*Fhatur0P(:,ixeq);
+ Fhat(:,ixeq) = Viconst{ixeq}*(gihat + HLV*(ci-Lit*gihat));
+ end
+ end
+ Fhat
+ Bhat = Fhat/A0hat; % ncoef-by-nvar reduced form lagged parameters.
+ end
+ nn = [nvar options_.ms.nlags imstp];
+ imfhat = fn_impulse(Bhat,swish,nn); % in the form that is congenial to RATS
+ imf3hat=reshape(imfhat,size(imfhat,1),nvar,nvar);
+ % imf3: row--steps, column--nvar responses, 3rd dimension--nvar shocks
+ imf3shat=permute(imf3hat,[1 3 2]);
+ % imf3s: permuted so that row--steps, column--nvar shocks,
+ % 3rd dimension--nvar responses
+ % Note: reshape(imf3s(1,:,:),nvar,nvar) = A0in (columns -- equations)
+ if options_.ms.indxgimfhat
+ figure
+ end
+ scaleout = fn_imcgraph(imfhat,nvar,imstp,xlab,ylab,options_.ms.indxgimfhat);
+ imfstd = max(abs(scaleout)'); % row: nvar (largest number); used for standard deviations
+
+ %
+ % %**** save stds. of both data and impulse responses in idfile1
+ % temp = [std(yactqmyge(:,3:end)); std(yactyrge(:,3:end)); imfstd]; %<<>>
+ % save idenyimstd.prn temp -ascii % export forecast and impulse response to the file "idenyimstd.prn", 3-by-nvar
+ % %
+ % %**** save stds. of both data and impulse responses in idfile1
+ % temp = [std(yactqmyge(:,3:end)); std(yactyrge(:,3:end)); imfstd]; %<<>>
+ % save idenyimstd.prn temp -ascii % export forecast and impulse response to the file "idenyimstd.prn", 3-by-nvar
+ % if options_.ms.indxparr
+ % idfile1='idenyimstd';
+ % end
+
+ %=====================================
+ % Now, out-of-sample forecasts. Note: Hm1t does not change with A0.
+ %=====================================
+ %
+ % * updating the last row of X (phi) with the current (last row of) y.
+ tcwx = nvar*options_.ms.nlags ; % total coefficeint without exogenous variables
+ phil = phi(size(phi,1),:);
+ phil(nvar+1:tcwx) = phil(1:tcwx-nvar);
+ phil(1:nvar) = y(end,:);
+ %*** exogenous variables excluding constant terms
+ if (nexo>1)
+ Xexoe = fn_dataext([yrEnd qmEnd],[yrEnd qmEnd],xdatae(:,[1:2 2+nvar+1:2+nvar+nexo-1]));
+ phil(1,tcwx+1:tcwx+nexo-1) = Xexoe(1,3:end);
+ end
+ %
+ %*** ML unconditional point forecast
+ nn = [nvar options_.ms.nlags nfqm];
+ if nexo<2
+ yforehat = fn_forecast(Bhat,phil,nn); % nfqm-by-nvar, in log
+ else
+ Xfexoe = fn_dataext(fdates(1,:),fdates(numel(fdates),:),xdatae(:,[1:2 2+nvar+1:2+nvar+nexo-1]));
+ %Xfexoe = fn_dataext(fdates(1,:),fdates(end,:),xdatae(:,[1:2 2+nvar+1:2+nvar+nexo-1]));
+ yforehat = fn_forecast(Bhat,phil,nn,nexo,Xfexoe(:,3:end)); % nfqm-by-nvar, in log
+ end
+ yforehate = [fdates yforehat];
+ %
+ yact1e = fn_dataext([yrEnd-nayr 1],[yrEnd qmEnd],xdatae(:,1:nvar+2));
+ if options_.ms.real_pseudo_forecast
+ %yact2e = fn_dataext([yrEnd-nayr 1],E2yrqm,xdatae);
+ yact2e = fn_dataext([yrEnd-nayr 1],[fdates(end,1) options_.ms.freq],xdatae(:,1:nvar+2));
+ else
+ yact2e=yact1e;
+ end
+ yafhate = [yact1e; yforehate]; % actual and forecast
+ %
+ %===== Converted to mg, qg, and calendar yg
+ %
+ [yafyrghate,yafyrhate,yafqmyghate] = fn_datana(yafhate,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno));
+ % actual and forecast growth rates
+ [yact2yrge,yact2yre,yact2qmyge] = fn_datana(yact2e,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno));
+ % only actual growth rates
+ yafyrghate
+ if options_.ms.indxgforhat
+ keyindx = [1:nvar];
+ conlab=['unconditional'];
+
+ figure
+ yafyrghate(:,3:end) = yafyrghate(:,3:end)/100;
+ yact2yrge(:,3:end) = yact2yrge(:,3:end)/100;
+ fn_foregraph(yafyrghate,yact2yrge,keyindx,rnum,cnum,options_.ms.freq,ylab,forelabel,conlab)
+ end
+
+ %-------------------------------------------------
+ % Setup for point conditional forecast
+ % ML Conditional Forecast
+ %-------------------------------------------------
+ %
+ %% See Zha's note "Forecast (1)" p. 5, RATS manual (some errors in RATS), etc.
+ %
+ %% Some notations: y(t+1) = y(t)B1 + e(t+1)inv(A0). e(t+1) is 1-by-n.
+ %% Let r(t+1)=e(t+1)inv(A0) + e(t+2)C + .... where inv(A0) is impulse
+ %% response at t=1, C at t=2, etc. The row of inv(A0) or C is
+ %% all responses to one shock.
+ %% Let r be q-by-1 (such as r(1) = r(t+1)
+ %% = y(t+1) (constrained) - y(t+1) (forecast)).
+ %% Use impulse responses to find out R (k-by-q) where k=nvar*nsteps
+ %% where nsteps the largest constrained step. The key of the program
+ %% is to creat R using impulse responses
+ %% Optimal solution for shock e where R'*e=r and e is k-by-1 is
+ %% e = R*inv(R'*R)*r.
+ %
+
+ if (nconstr > 0)
+ %*** initializing
+ stepcon=cell(nconstr,1); % initializing, value y conditioned
+ valuecon=zeros(nconstr,1); % initializing, value y conditioned
+ varcon=zeros(nconstr,1); % initializing, endogous variables conditioned
+ varcon(:)=PorR; % 1: Pcm; 3: FFR; 5: CPI
+
+ %
+ for i=1:nconstr
+ if i<=nconstr1
+ stepcon{i}=i; % FFR
+ elseif i<=2*nconstr1
+ stepcon{i}=i-nconstr1; % FFR
+ elseif i<=3*nconstr1
+ stepcon{i}=i-2*nconstr1; % FFR
+ elseif i<=4*nconstr1
+ stepcon{i}=i-3*nconstr1; % FFR
+ elseif i<=5*nconstr1
+ stepcon{i}=i-4*nconstr1; % FFR
+ elseif i<=6*nconstr1
+ stepcon{i}=i-5*nconstr1; % FFR
+ end
+ end
+
+% for i=1:nconstr
+% stepcon{i}=i; % FFR
+% end
+
+% bend=12;
+% stepcon{1}=[1:bend]'; % average over
+% stepcon{nconstr1+1}=[1:options_.ms.freq-qmSub]'; % average over the remaing months in 1st forecast year
+% stepcon{nconstr1+2}=[options_.ms.freq-qmSub+1:options_.ms.freq-qmSub+12]'; % average over 12 months next year
+% stepcon{nconstr1+3}=[options_.ms.freq-qmSub+13:options_.ms.freq-qmSub+24]'; % average over 12 months. 3rd year
+% stepcon{nconstr1+4}=[options_.ms.freq-qmSub+25:options_.ms.freq-qmSub+36]'; % average over 12 months. 4th year
+
+% %**** avearage condition over, say, options_.ms.freq periods
+% if qmEnd==options_.ms.freq
+% stepcon{1}=[1:options_.ms.freq]'; % average over the remaing periods in 1st forecast year
+% else
+% stepcon{1}=[1:options_.ms.freq-qmEnd]'; % average over the remaing periods in 1st forecast year
+% end
+% for kj=2:nconstr
+% stepcon{kj}=[length(stepcon{kj-1})+1:length(stepcon{kj-1})+options_.ms.freq]'; % average over 12 months next year
+% end
+
+ if options_.ms.real_pseudo_forecast
+% %*** conditions in every period
+% for i=1:nconstr
+% valuecon(i) = yact(actup+i,varcon(i));
+% %valuecon(i) = mean( yact(actup+1:actup+bend,varcon(i)) );
+% %valuecon(i) = 0.060; % 95:01
+% %valuecon(i) = (0.0475+0.055)/2; % 94:10
+% end
+
+% %*** average condtions over,say, options_.ms.freq periods.
+% for i=nconstr1+1:nconstr1+nconstr2
+% i=1;
+% valuecon(nconstr1+i) = ( ( mean(ylast12Cal(:,varcon(nconstr1+i)),1) + ...
+% log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100) )*options_.ms.freq - ...
+% yCal_1(:,varcon(nconstr1+i)) ) ./ length(stepcon{nconstr1+i});
+% % the same as unconditional "yactCalyg" 1st calendar year
+% i=2;
+% valuecon(nconstr1+i) = mean(ylast12Cal(:,varcon(nconstr1+i))) + ...
+% log(1+yactCalyg(yAg-yFg+1,varcon(nconstr1+i))/100) ...
+% + log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100);
+% % the same as actual "yactCalgy" 2nd calendar year
+% i=3;
+% valuecon(nconstr1+i) = valuecon(nconstr1+i-1) + ...
+% log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100);
+% % the same as actual "yactCalgy" 3rd calendar year
+% %i=4;
+% %valuecon(nconstr1+i) = valuecon(nconstr1+i-1) + ...
+% % log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100);
+% % the same as actual "yactCalgy" 4th calendar year
+% end
+
+ %*** conditions in every period
+ vpntM = fn_dataext(E1yrqm, E2yrqm,xdatae); % point value matrix with dates
+ % vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates
+ for i=1:nconstr
+ if i<=nconstr1
+ valuecon(i) = vpntM(i,2+varcon(i)); % 2: first 2 elements are dates
+ elseif i<=2*nconstr1
+ valuecon(i) = vpntM(i-nconstr1,2+varcon(i));
+ elseif i<=3*nconstr1
+ valuecon(i) = vpntM(i-2*nconstr1,2+varcon(i));
+ elseif i<=4*nconstr1
+ valuecon(i) = vpntM(i-3*nconstr1,2+varcon(i));
+ elseif i<=5*nconstr1
+ valuecon(i) = vpntM(i-4*nconstr1,2+varcon(i));
+ elseif i<=6*nconstr1
+ valuecon(i) = vpntM(i-5*nconstr1,2+varcon(i));
+ end
+ end
+
+% %*** average condtions over,say, options_.ms.freq periods.
+% if qmEnd==options_.ms.freq
+% vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates
+% valuecon(1) = vaveM(1,2+varcon(1)); % 2: first 2 elements are dates
+% else
+% vaveM = fn_dataext([yrEnd 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates
+% yactrem = fn_dataext([yrEnd qmEnd+1],[yrEnd options_.ms.freq],xdatae);
+% valuecon(1) = sum(yactrem(:,2+varcon(1)),1)/length(stepcon{1});
+% % 2: first 2 elements are dates
+% end
+% for kj=2:nconstr
+% valuecon(kj) = vaveM(kj,2+varcon(kj)); % 2: first 2 elements are dates
+% end
+ else
+ vpntM = dataext([yrEnd qmEnd+1],[yrEnd qmEnd+2],xdatae); % point value matrix with dates
+ for i=1:nconstr
+ if i<=nconstr1
+ valuecon(i) = vpntM(i,2+varcon(i)); % 2: first 2 elements are dates; Poil
+ elseif i<=2*nconstr1
+ valuecon(i) = vpntM(i-nconstr1,2+varcon(i)); % 2: first 2 elements are dates; M2
+ elseif i<=3*nconstr1
+ valuecon(i) = vpntM(i-2*nconstr1,2+varcon(i)); % 2: first 2 elements are dates; FFR
+ elseif i<=4*nconstr1
+ valuecon(i) = vpntM(i-3*nconstr1,2+varcon(i)); % 2: first 2 elements are dates; CPI
+ elseif i<=5*nconstr1
+ valuecon(i) = vpntM(i-4*nconstr1,2+varcon(i)); % 2: first 2 elements are dates; U
+ elseif i<=5*nconstr1+nconstr2
+ valuecon(i)=xdata(end,5)+(i-5*nconstr1)*log(1.001)/options_.ms.freq; %CPI
+ elseif i<=5*nconstr1+2*nconstr2
+ valuecon(i)=0.0725; %FFR
+ else
+ valuecon(i)=xdata(end,6)+(i-5*nconstr1-2*nconstr2)*0.01/nfqm; %U
+ end
+ end
+ %valuecon(i) = 0.060; % 95:01
+ end
+ else
+ valuecon = [];
+ stepcon = [];
+ varcon = [];
+ end
+
+ nstepsm = 0; % initializing, the maximum step in all constraints
+ for i=1:nconstr
+ nstepsm = max([nstepsm max(stepcon{i})]);
+ end
+
+ if cnum
+ if options_.ms.real_pseudo_forecast && options_.ms.banact
+ for i=1:length(banstp{1})
+ banval{1}(1:length(banstp{1}),1) = ...
+ yactCalyg(yAg-yFg+1:yAg-yFg+length(banstp{1}),banvar(1)) - 2;
+ banval{1}(1:length(banstp{1}),2) = ...
+ yactCalyg(yAg-yFg+1:yAg-yFg+length(banstp{1}),banvar(1)) + 2;
+ end
+ end
+ end
+
+
+ %===================================================
+ % ML conditional forecast
+ %===================================================
+ %/*
+ [ychat,Estr,rcon] = fn_fcstidcnd(valuecon,stepcon,varcon,nstepsm,...
+ nconstr,options_.ms.eq_ms,nvar,options_.ms.nlags ,phil,0,0,yforehat,imf3shat,A0hat,Bhat,...
+ nfqm,options_.ms.tlindx,options_.ms.tlnumber,options_.ms.ncms,options_.ms.eq_cms);
+ ychate = [fdates ychat];
+ yachate = [yact1e; ychate]; % actual and condtional forecast
+ %===== Converted to mg, qg, and calendar yg
+ [yacyrghate,yacyrhate,yacqmyghate] = fn_datana(yachate,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno));
+ % actual and conditional forecast growth rates
+ if options_.ms.indxgdls && nconstr
+ keyindx = [1:nvar];
+ % conlab=['conditional on' ylab{PorR(1)}];
+ conlab=['v-conditions'];
+
+ figure
+ fn_foregraph(yafyrghate,yact2yrge,keyindx,rnum,cnum,options_.ms.freq,ylab,forelabel,conlab)
+ end
+
+ if options_.ms.ncsk
+ Estr = zeros(nfqm,nvar);
+ Estr(1:2,:) = [
+ -2.1838 -1.5779 0.53064 -0.099425 -0.69269 -1.0391
+ 1.9407 3.3138 -0.10563 -0.55457 -0.68772 1.3534
+ ];
+ Estr(3:6,3) = [0.5*ones(1,4)]'; % MD shocks
+
+ Estr(3:10,2) = [1.5 1.5 1.5*ones(1,6)]'; % MS shocks
+
+ %Estr(3:6,6) = 1*ones(4,1); % U shocks
+ %Estr(8:11,4) = 1*ones(4,1); % y shocks
+
+ %Estr(3:10,2) = [2.5 2.5 1.5*ones(1,6)]'; % MS shocks alone
+
+ nn = [nvar options_.ms.noptions_.ms.nlags nfqm];
+ ycEhat = forefixe(A0hat,Bhat,phil,nn,Estr);
+ ycEhate = [fdates ycEhat];
+ yacEhate = [yact1e; ycEhate]; % actual and condtional forecast
+ %===== Converted to mg, qg, and calendar yg
+ [yacEyrghate,yacEyrhate,yacEqmyghate] = datana(yacEhate,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno));
+ % actual and conditional forecast growth rates
+ disp([sprintf('%4.0f %2.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',yacEyrghate')])
+
+ if 1
+ keyindx = [1:nvar];
+ % conlab=['conditional on' ylab{PorR(1)}];
+ conlab=['shock-conditions'];
+
+ figure
+ gyrfore(yacEyrghate,yact2yrge,keyindx,rnum,cnum,ylab,forelabel,conlab)
+ end
+ end
+
+ %-----------------------------------------------------------
+ % Compute structural shocks for the whole sample period excluding dummy observations.
+ %-----------------------------------------------------------
+ ywod = y(options_.ms.dummy_obs +1:end,:); % without dummy observations
+ phiwod=phi(options_.ms.dummy_obs +1:end,:); % without dummy observations
+ eplhat=ywod*A0hat-phiwod*Fhat;
+ qmStartWod = mod(qmStart+options_.ms.nlags ,options_.ms.freq);
+ if (~qmStartWod)
+ qmStartWod = options_.ms.freq;
+ end
+ yrStartWod = yrStart + floor((qmStart+options_.ms.nlags -1)/options_.ms.freq);
+ dateswod = fn_dataext([yrStartWod qmStartWod],[yrEnd qmEnd],xdatae(:,[1:2]));
+ eplhate = [dateswod eplhat];
+
+ Aphat = Fhat;
+end
+
+%----------------------------------------
+% Tests for LR, HQ, Akaike, Schwarz. The following gives a guidance.
+% But the computation has to be done in a different M file by exporting fhat's
+% from different idfile's.
+%----------------------------------------
+%
+%if ~options_.ms.contemp_reduced_form
+% SpHR=A0in'*A0in;
+%end
+%%
+%if ~isnan(SpHR) && ~options_.ms.contemp_reduced_form
+% warning(' ')
+% disp('Make sure you run the program with options_.ms.contemp_reduced_form =1 first.')
+% disp('Otherwise, the following test results such as Schwartz are incorrect.')
+% disp('All other results such as A0ml and imfs, however, are correct.')
+% disp('Press anykey to contintue or ctrl-c to stop now')
+% pause
+
+% load SpHUout
+
+% logLHU=-fss*sum(log(diag(chol(SpHU)))) -0.5*fss*nvar % unrestricted logLH
+
+% logLHR=-fhat % restricted logLH
+% tra = reshape(SpHU,nvar*nvar,1)'*reshape(A0*A0',nvar*nvar,1);
+% df=(nvar*(nvar+1)/2 - length(a0indx));
+% S=2*(logLHU-logLHR);
+% SC = (nvar*(nvar+1)/2 - length(a0indx)) * log(fss);
+% disp(['T -- effective sample size: ' num2str(fss)])
+% disp(['Trace in the overidentified posterior: ' num2str(tra)])
+% disp(['Chi2 term -- 2*(logLHU-logLHR): ' num2str(S)])
+% disp(['Degrees of freedom: ' num2str(df)])
+% disp(['SC -- df*log(T): ' num2str(SC)])
+% disp(['Akaike -- 2*df: ' num2str(2*df)])
+% disp(['Classical Asymptotic Prob at chi2 term: ' num2str(cdf('chi2',S,df))])
+
+% %*** The following is the eigenanalysis in the difference between
+% %*** unrestricted (U) and restricted (R)
+% norm(A0'*SpHU*A0-diag(diag(ones(6))))/6;
+% norm(SpHU-A0in'*A0in)/6;
+
+% corU = corr(SpHU);
+% corR = corr(SpHR);
+
+% [vU,dU]=eigsort(SpHU,1);
+% [vR,dR]=eigsort(SpHR,1);
+
+% [log(diag(dU)) log(diag(dR)) log(diag(dU))-log(diag(dR))];
+
+% sum(log(diag(dU)));
+% sum(log(diag(dR)));
+%else
+% disp('To run SC test, turn options_.ms.contemp_reduced_form =1 first and then turn options_.ms.contemp_reduced_form =0')
+%end
+
+
+%***** Simply regression
+%X=[phi(:,3) y(:,2)-phi(:,2) y(:,1)-phi(:,7) ones(fss,1)];
+%� Y=y(:,3);
+%� b=regress(Y,X)
+
+%=== Computes the roots for the whole system.
+rootsinv = fn_varoots(Bhat,nvar,options_.ms.nlags )
+abs(rootsinv)
+
+
+bhat =xhat;
+n0const=n0; % For constant parameter models.
+n0const=n0; % For constant parameter models.
+npconst=np; % For constant parameter models.
+save outdata_a0dp_const.mat A0hat bhat Aphat n0const
diff --git a/matlab/ms-sbvar/msstart_setup.m b/matlab/ms-sbvar/msstart_setup.m
index 31f10ca90..3430d98c2 100644
--- a/matlab/ms-sbvar/msstart_setup.m
+++ b/matlab/ms-sbvar/msstart_setup.m
@@ -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 .
-
-% ** 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> 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 .
+
+% ** 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> 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
+%<<<<<<<<<<<<<<<<<<<
+
+
+
+
+