ms_sbvar_setup

 function ms_sbvar_setup(options_)
 does the general file initialization for ms sbvar

 INPUTS
    options_:    (struct)    options

 OUTPUTS
    none

 SPECIAL REQUIREMENTS
    none

0001 function ms_sbvar_setup(options_)
0002 % function ms_sbvar_setup(options_)
0003 % does the general file initialization for ms sbvar
0004 %
0005 % INPUTS
0006 %    options_:    (struct)    options
0007 %
0008 % OUTPUTS
0009 %    none
0010 %
0011 % SPECIAL REQUIREMENTS
0012 %    none
0013 
0014 % Copyright (C) 2003-2011 Dynare Team
0015 %
0016 % This file is part of Dynare.
0017 %
0018 % Dynare is free software: you can redistribute it and/or modify
0019 % it under the terms of the GNU General Public License as published by
0020 % the Free Software Foundation, either version 3 of the License, or
0021 % (at your option) any later version.
0022 %
0023 % Dynare is distributed in the hope that it will be useful,
0024 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0025 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0026 % GNU General Public License for more details.
0027 %
0028 % You should have received a copy of the GNU General Public License
0029 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0030 
0031 if ~isfield(options_.ms,'initial_year')
0032     error('Must set initial_year option');
0033 end
0034 
0035 if ~isfield(options_.ms,'final_year')
0036     error('Must set final_year option');
0037 end
0038 
0039 if ~isfield(options_,'datafile')
0040     error('Must set datafile option');
0041 end
0042 
0043 options_.data = read_variables(options_.datafile, ...
0044     options_.varobs, [], options_.xls_sheet, options_.xls_range);
0045 
0046 if options_.ms.upper_cholesky
0047     if options_.ms.lower_cholesky
0048         error(['Upper Cholesky and lower Cholesky decomposition can''t be ' ...
0049                'requested at the same time!'])
0050     else
0051         options_.ms.restriction_fname = 'upper_cholesky';
0052     end
0053 elseif options_.ms.lower_cholesky
0054     options_.ms.restriction_fname = 'lower_cholesky';
0055 elseif ~isempty(options_.ms.Qi) && ~isempty(options_.ms.Ri)
0056     options_.ms.restriction_fname = 'exclusions';
0057 else
0058     options_.ms.restriction_fname = 0;
0059 end
0060 
0061 %==========================================================================
0062 %== Markov Process Specification File
0063 %==========================================================================
0064 markov_file = [options_.ms.output_file_tag '_markov_file.dat'];
0065 
0066 %==========================================================================
0067 %== BVAR prior
0068 %==========================================================================
0069 
0070 %=== The following mu is effective only if indxPrior==1.
0071 %mu = zeros(6,1);   % hyperparameters
0072 if length(options_.ms.coefficients_prior_hyperparameters) ~= 6
0073     error('When specifying the coefficients_prior_hyperparameters, you must pass a vector of 6 numbers')
0074 end
0075 mu = options_.ms.coefficients_prior_hyperparameters;
0076 mu = reshape(mu,1,6);
0077 
0078 % mu(1): overall tightness for A0 and Aplus
0079 % mu(2): relative tightness for Aplus
0080 % mu(3): relative tightness for the constant term
0081 % mu(4): tightness on lag decay.  (1.2 - 1.5 faster decay produces better
0082 % inflation forecasts
0083 % mu(5): weight on nvar sums of coeffs dummy observations (unit roots).
0084 % mu(6): weight on single dummy initial observation including constant
0085 % (cointegration, unit roots, and stationarity).
0086 
0087 % Alpha on p. 66 for squared time-varying structural shock lambda.
0088 galp = options_.ms.alpha;
0089 
0090 % Beta on p. 66 for squared time-varying structural shock lambda.
0091 gbeta = options_.ms.beta;
0092 
0093 % Case 3 (no state change across options_.ms.nlags (l) but allows all variables for a
0094 % given lag to switch states). Normal prior variance for glamda
0095 % (nvar-by-nvar for each state) for different variables in lagged D+.  See
0096 % p.71v.
0097 gsig2_lmdm = options_.ms.gsig2_lmdm;
0098 
0099 
0100 %==========================================================================
0101 %== Data
0102 %==========================================================================
0103 % Read in data to produce rectangular array named xdd.  Each column is one
0104 % data series.
0105 xdd=options_.data;
0106 
0107 % Information about timing of the data for consistancy checks
0108 % quarters (4) or months (12)
0109 q_m = options_.ms.freq;
0110 % beginning year in data set
0111 yrBin=options_.ms.initial_year;
0112 % beginning quarter or month in data set
0113 %options_.ms.initial_subperiod = 1;
0114 qmBin=options_.ms.initial_subperiod;
0115 % final year in data set
0116 yrFin=options_.ms.final_year;
0117 % final month or quarter in data set
0118 qmFin=options_.ms.final_subperiod;
0119 % first year to use in estimation
0120 yrStart=options_.ms.initial_year;
0121 % first quarter or month to use in estimation
0122 qmStart=options_.ms.initial_subperiod;
0123 % last year to use in estimation
0124 yrEnd=options_.ms.final_year;
0125 % last quater or month to use in estimation
0126 qmEnd=options_.ms.final_subperiod;
0127 % Log variables in xdd
0128 logindx = [];
0129 
0130 % Convert percent to decimal in xdd
0131 pctindx = [];
0132 
0133 % Select the variable to use and rearrange columns if desired
0134 %vlist = [3 1 2];
0135 %options_.ms.vlist = [1 2 3];
0136 options_.ms.vlist = 1:size(options_.varobs,1);
0137 vlist1=options_.ms.vlist;
0138 
0139 %==========================================================================
0140 %==========================================================================
0141 %==========================================================================
0142 %== Beginning of code.  Modify below at own risk.
0143 %==========================================================================
0144 
0145 % options that may at some point become user specified
0146 %indxC0Pres = 0;    % 1: cross-A0-and-A+ restrictions; 0: idfile_const is all we have
0147 indxC0Pres =options_.ms.cross_restrictions;
0148 % Example for indxOres==1: restrictions of the form P(t) = P(t-1).
0149 %Rform = 0;         % 1: contemporaneous recursive reduced form; 0: restricted (non-recursive) form
0150 Rform =options_.ms.contemp_reduced_form;
0151 % % % Pseudo = 0;  % 1: Pseudo forecasts; 0: real time forecasts
0152 %indxPrior = 1;     % 1: Bayesian prior; 0: no prior
0153 indxPrior =options_.ms.bayesian_prior;
0154 %indxDummy = indxPrior;  % 1: add dummy observations to the data; 0: no dummy added.
0155 indxDummy = options_.ms.bayesian_prior;
0156 %ndobs = 0;         % No dummy observations for xtx, phi, fss, xdatae, etc.  Dummy observations are used as an explicit prior in fn_rnrprior_covres_dobs.m.
0157 %ndobs =options_.ms.dummy_obs;
0158 %if indxDummy
0159 %   ndobs=nvar+1;         % number of dummy observations
0160 %else
0161 %   ndobs=0;    % no dummy observations
0162 %end
0163 
0164 %
0165 hpmsmd = [0.0; 0.0];
0166 indxmsmdeqn = [0; 0; 0; 0];  %This option disenable using this in fn_rnrprior_covres_dobs.m
0167 
0168 nStates = -1;
0169 
0170 %==========================================================================
0171 %== Create initialization file
0172 %==========================================================================
0173 
0174 %======================================================================
0175 %== Check and setup data
0176 %======================================================================
0177 % log data
0178 xdd(:,logindx) = log(xdd(:,logindx));
0179 
0180 % convert percentage to decimal
0181 xdd(:,pctindx)=.01*xdd(:,pctindx);
0182 
0183 if (q_m ~= 12) && (q_m ~= 4)
0184     disp('Warning: data must be monthly or quarterly!')
0185     return
0186 end
0187 
0188 % number of data points
0189 nData=(yrFin-yrBin)*q_m + (qmFin-qmBin+1);
0190 % number of data points in estimation sample
0191 nSample=(yrEnd-yrStart)*q_m + (qmEnd-qmEnd+1);
0192 % number of periods not used at beginning of data (non-negative number)
0193 nStart=(yrStart-yrBin)*q_m + (qmStart-qmBin);
0194 % number of periods not used at end of data (non-positive number)
0195 nEnd=(yrEnd-yrFin)*q_m + (qmEnd-qmFin);
0196 
0197 if (nEnd > 0) || (nStart < 0)
0198     disp('Warning: desired estimation period not in data set!')
0199     return
0200 end
0201 if (nSample <= 0)
0202     disp('Warning: no data points in estimation period!')
0203     return
0204 end
0205 
0206 % reorder variables and create estimation data set
0207 xdgel=xdd(nStart+1:nData+nEnd,vlist1);
0208 
0209 % bad data points
0210 baddata = find(isnan(xdgel));
0211 if ~isempty(baddata)
0212     disp('Warning: some data for estimation period are unavailable.')
0213     return
0214 end
0215 
0216 % set nvar and nexo
0217 nvar=size(xdgel,2);
0218 nexo=1;
0219 
0220 % Arranged data information, WITHOUT dummy obs when 0 after mu is used.
0221 % See fn_rnrprior_covres_dobs.m for using the dummy observations as part of
0222 % an explicit prior.
0223 [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh] = fn_dataxy(nvar,options_.ms.nlags,xdgel,mu,0,nexo);
0224 
0225 
0226 %======================================================================
0227 %== Linear Restrictions
0228 %======================================================================
0229 if Rform
0230     Ui=cell(nvar,1); Vi=cell(ncoef,1);
0231     for kj=1:nvar
0232         Ui{kj} = eye(nvar);  Vi{kj} = eye(ncoef);
0233     end
0234 else
0235     [Ui,Vi,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,nvar,nexo,options_.ms);
0236     if min(n0)==0
0237         disp('A0: restrictions give no free parameters in one of equations')
0238         return
0239     elseif min(np)==0
0240         disp('Aplus: Restrictions in give no free parameters in one of equations')
0241         return
0242     end
0243 end
0244 
0245 
0246 %======================================================================
0247 %== Estimation
0248 %======================================================================
0249 if indxPrior
0250     %*** Obtains asymmetric prior (with no linear restrictions) with dummy observations as part of an explicit prior (i.e,
0251     %      reflected in Hpmulti and Hpinvmulti).  See Forecast II, pp.69a-69b for details.
0252     if 1  % Liquidity effect prior on both MS and MD equations.
0253         [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior_covres_dobs(nvar,q_m,options_.ms.nlags,xdgel,mu,indxDummy,hpmsmd,indxmsmdeqn);
0254     else
0255         [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior(nvar,q_m,options_.ms.nlags,xdgel,mu);
0256     end
0257     
0258     %*** Combines asymmetric prior with linear restrictions
0259     [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Ui,Vi,Pi,H0multi,Hpmulti,nvar);
0260     
0261     %*** Obtains the posterior matrices for estimation and inference
0262     [Pmat,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0invtld,Hpinvtld,Ui,Vi);
0263 else
0264     %*** Obtain the posterior matrices for estimation and inference
0265     [Pmat,H0inv,Hpinv] = fn_dlrpostr(xtx,xty,yty,Ui,Vi);
0266 end
0267 
0268 if Rform
0269     %*** Obtain the ML estimate
0270     A0hatinv = chol(H0inv{1}/fss);   % upper triangular but lower triangular choleski
0271     A0hat=inv(A0hatinv);
0272     
0273     Aphat = Pmat{1}*A0hat;
0274 else
0275     %*** Obtain the ML estimate
0276     %   load idenml
0277     x = 10*rand(sum(n0),1);
0278     H0 = eye(sum(n0));
0279     crit = 1.0e-9;
0280     nit = 10000;
0281     %
0282     [fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Ui,nvar,n0,fss,H0inv);
0283 
0284     A0hat = fn_tran_b2a(xhat,Ui,nvar,n0);
0285     
0286     xhat = fn_tran_a2b(A0hat,Ui,nvar,n0);
0287     [Aphat,ghat] = fn_gfmean(xhat,Pmat,Vi,nvar,ncoef,n0,np);
0288     if indxC0Pres
0289         Fhatur0P = Fhat;  % ur: unrestriced across A0 and A+
0290         for ki = 1:size(ixmC0Pres,1)   % loop through the number of equations in which
0291             % cross-A0-A+ restrictions occur. See St. Louis Note p.5.
0292             ixeq = ixmC0Pres{ki}(1,1);   % index for the jth equation in consideration.
0293             Lit = Vi{ixeq}(ixmC0Pres{ki}(:,2),:);  % transposed restriction matrix Li
0294             % V_j(i,:) in f_j(i) = V_j(i,:)*g_j
0295             ci = ixmC0Pres{ki}(:,4) .* A0hat(ixmC0Pres{ki}(:,3),ixeq);
0296             % s * a_j(h) in the restriction f_j(i) = s * a_j(h).
0297             LtH = Lit/Hpinv{ixeq};
0298             HLV = LtH'/(LtH*Lit');
0299             gihat = Vi{ixeq}'*Fhatur0P(:,ixeq);
0300             Aphat(:,ixeq) = Vi{ixeq}*(gihat + HLV*(ci-Lit*gihat));
0301         end
0302     end
0303 end
0304 
0305 
0306 %======================================================================
0307 %== Create matlab initialization file
0308 %======================================================================
0309 matlab_filename = ['matlab_',options_.ms.output_file_tag,'.prn'];
0310 fidForC = fopen(matlab_filename,'w');
0311 
0312 fprintf(fidForC,'\n%s\n','//== gxia: alpha parameter for gamma prior of xi ==//');
0313 fprintf(fidForC,' %20.15f ', galp);
0314 fprintf(fidForC, '\n\n');
0315 
0316 fprintf(fidForC,'\n%s\n','//== gxib: beta parameter for gamma prior of xi ==//');
0317 fprintf(fidForC,' %20.15f ', gbeta);
0318 fprintf(fidForC, '\n\n');
0319 
0320 fprintf(fidForC,'\n%s\n','//== glamdasig: sigma parameter for normal prior of lamda ==//');
0321 fprintf(fidForC,' %20.15f ', sqrt(gsig2_lmdm));
0322 fprintf(fidForC, '\n\n');
0323 
0324 %=== lags, nvar, nStates, sample size (excluding options_.ms.nlags where, with dummyies, fss=nSample-options_.ms.nlags+ndobs).
0325 fprintf(fidForC,'\n%s\n','//== lags, nvar, nStates, T ==//');
0326 fprintf(fidForC,' %d  %d  %d  %d\n\n\n',options_.ms.nlags, nvar, nStates, fss);
0327 
0328 %=== A0hat nvar-by-nvar from the constant VAR.
0329 fprintf(fidForC,'\n%s\n','//== A0hat: nvar-by-nvar ==//');
0330 indxFloat = 1;
0331 xM = A0hat;
0332 nrows = nvar;
0333 ncols = nvar;
0334 fn_fprintmatrix(fidForC, xM, nrows, ncols, indxFloat)
0335 
0336 %=== Aphat ncoef-by-nvar from the constant VAR.
0337 %=== Each column of Aphat is in the order of [nvar variables for 1st lag, ..., nvar variables for last lag, constant term].
0338 fprintf(fidForC,'\n%s\n','//== Aphat: ncoef(lags*nvar+1)-by-nvar ==//');
0339 indxFloat = 1;
0340 xM = Aphat;
0341 nrows = ncoef;
0342 ncols = nvar;
0343 fn_fprintmatrix(fidForC, xM, nrows, ncols, indxFloat)
0344 
0345 %=== n0const: nvar-by-1, whose ith element represents the number of free A0 parameters in ith equation for the case of constant parameters.
0346 fprintf(fidForC,'\n%s\n','//== n0const: nvar-by-1 ==//');
0347 indxFloat = 0;
0348 xM = n0;
0349 nrows = 1;
0350 ncols = nvar;
0351 fn_fprintmatrix(fidForC, xM', nrows, ncols, indxFloat)
0352 
0353 %=== npconst: nvar-by-1, whose ith element represents the number of free A+ parameters in ith equation for the case of constant parameters.
0354 fprintf(fidForC,'\n%s\n','//== npconst: nvar-by-1 ==//');
0355 indxFloat = 0;
0356 xM = np;
0357 nrows = 1;
0358 ncols = nvar;
0359 fn_fprintmatrix(fidForC, xM', nrows, ncols, indxFloat)
0360 
0361 %=== Specification
0362 fprintf(fidForC,'\n%s','//== Specification (0=default  1=Sims-Zha  2=Random Walk) ==//');
0363 fprintf(fidForC,'\n%d\n\n',options_.ms.specification);
0364 
0365 %=== Uiconst: nvar-by-1 cell.  In each cell, nvar-by-qi orthonormal basis for the null of the ith
0366 %           equation contemporaneous restriction matrix where qi is the number of free parameters.
0367 %           With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
0368 %           of total original parameters and bi is a vector of free parameters. When no
0369 %           restrictions are imposed, we have Ui = I.  There must be at least one free
0370 %           parameter left for the ith equation.
0371 fprintf(fidForC,'\n%s\n','//== Uiconst: cell(nvar,1) and nvar-by-n0const(i) for the ith cell (equation) ==//');
0372 for i_=1:nvar
0373     fn_fprintmatrix(fidForC, Ui{i_}, nvar, n0(i_), 1);
0374 end
0375 
0376 %=== Viconst: nvar-by-1 cell.  In each cell, k-by-ri orthonormal basis for the null of the ith
0377 %           equation lagged restriction matrix where k is a total of exogenous variables and
0378 %           ri is the number of free parameters. With this transformation, we have fi = Vi*gi
0379 %           or Vi'*fi = gi where fi is a vector of total original parameters and gi is a
0380 %           vector of free parameters. There must be at least one free parameter left for
0381 %           the ith equation.
0382 fprintf(fidForC,'\n%s\n','//== Viconst: cell(nvar,1) and ncoef-by-n0const(i) for the ith cell (equation) ==//');
0383 for i_=1:nvar
0384     fn_fprintmatrix(fidForC, Vi{i_}, ncoef, np(i_), 1);
0385 end
0386 
0387 %=== H0barconstcell: cell(nvar,1) (equations) and n-by-n for each cell (equaiton).
0388 %=== H0barconst:  prior covariance matrix for each column of A0 under asymmetric prior (including SZ dummy obs.) with NO linear restrictions imposed yet.
0389 fprintf(fidForC,'\n%s\n','//== H0barconstcell: cell(nvar,1) and n-by-n for the ith cell (equation) ==//');
0390 for i_=1:nvar
0391     fn_fprintmatrix(fidForC, H0multi(:,:,i_), nvar, nvar, 1);
0392 end
0393 
0394 %=== Hpbarconstcell: cell(nvar,1) (equations) and ncoef-by-ncoef for each cell (equaiton).
0395 %=== Hpbarconst:  prior covariance matrix for each column of A+ under asymmetric prior (including SZ dummy obs.) with NO linear restrictions imposed yet.
0396 fprintf(fidForC,'\n%s\n','//== Hpbarconstcell: cell(nvar,1) and ncoef-by-ncoef for the ith cell (equation) ==//');
0397 for i_=1:nvar
0398     fn_fprintmatrix(fidForC, Hpmulti(:,:,i_), ncoef, ncoef, 1);
0399 end
0400 
0401 %=== phi:  X; T-by-k; column: [nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term]
0402 fprintf(fidForC,'\n%s\n','//== Xright -- X: T-by-ncoef ==//');
0403 xM = phi;
0404 nrows = fss;
0405 ncols = ncoef;
0406 for ki=1:nrows
0407     for kj=1:ncols
0408         fprintf(fidForC,' %20.15f ',xM((kj-1)*nrows+ki));
0409         if (kj==ncols)
0410             fprintf(fidForC,'\n');
0411         end
0412     end
0413     if (ki==nrows)
0414         fprintf(fidForC,'\n\n');
0415     end
0416 end
0417 
0418 %=== y:    Y: T-by-nvar where T=fss
0419 fprintf(fidForC,'\n%s\n','//== Yleft -- Y: T-by-nvar ==//');
0420 xM = y;
0421 nrows = fss;
0422 ncols = nvar;
0423 for ki=1:nrows
0424     for kj=1:ncols
0425         fprintf(fidForC,' %20.15f ',xM((kj-1)*nrows+ki));
0426         if (kj==ncols)
0427             fprintf(fidForC,'\n');
0428         end
0429     end
0430     if (ki==nrows)
0431         fprintf(fidForC,'\n\n');
0432     end
0433 end
0434 
0435 fclose(fidForC);
0436 
0437 %======================================================================
0438 %== Create C initialization filename
0439 %======================================================================
0440 ms_write_markov_file(markov_file,options_)
0441 create_init_file = [matlab_filename,' ',markov_file,' ',options_.ms.file_tag];
0442 [err] = ms_sbvar_create_init_file(create_init_file);
0443 mexErrCheck('ms_sbvar_create_init_file',err);
0444 end