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 20040001 %function msstart2(options_) 0002 % This .m file is called for point graphics or error bands and 0003 % It starts for both data and Bayesian estimation (when IxEstimat==0, 0004 % no estimation but only data analysis), which allows you to set 0005 % (a) available data range, 0006 % (b) sample range, 0007 % (c) rearrangement of actual data such as mlog, qg, yg 0008 % (d) conditions of shocks 'Cms', 0009 % (c) conditions of variables 'nconstr', 0010 % (e) soft conditions 'nbancon,' 0011 % (f) produce point conditional forecast (at least conditions on variables). 0012 % 0013 % February 2004 0014 0015 % Copyright (C) 2011 Dynare Team 0016 % 0017 % This file is part of Dynare. 0018 % 0019 % Dynare is free software: you can redistribute it and/or modify 0020 % it under the terms of the GNU General Public License as published by 0021 % the Free Software Foundation, either version 3 of the License, or 0022 % (at your option) any later version. 0023 % 0024 % Dynare is distributed in the hope that it will be useful, 0025 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0026 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0027 % GNU General Public License for more details. 0028 % 0029 % You should have received a copy of the GNU General Public License 0030 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0031 0032 % ** ONLY UNDER UNIX SYSTEM 0033 %path(path,'/usr2/f1taz14/mymatlab') 0034 0035 %addpath('C:\SoftWDisk\MATLAB6p5\toolbox\cstz') 0036 0037 0038 msstart_setup 0039 %=========================================== 0040 % Exordium II 0041 %=========================================== 0042 %options_.ms.ncsk = 0; % conditional directly on shoocks. Unlike Cms, not on variables that back out shocks 0043 %options_.ms.nstd = 6; % number of standard deviations to cover the range in which distributions are put to bins 0044 %options_.ms.ninv = 1000; % the number of bins for grouping, say, impulse responses 0045 %options_.ms.Indxcol = [1:nvar]; % a vector of random columns in which MC draws are made. 0046 % 0047 %options_.ms.indxparr = 1; % 1: parameters random; 0: no randomness in parameters 0048 % Note, when 0, there is no effect from the values of options_.ms.IndxAp, options_.ms.Aband, etc. 0049 %options_.ms.indxovr = 0; % 1: distributions for other variables of interest; 0: no distribution. 0050 % Example: joint distribution of a(1) and a(2). Only for specific purposes 0051 %options_.ms.Aband = 1; % 1: error bands with only A0 and A+ random. 0052 %options_.ms.IndxAp = 1; % 1: generate draws of A+; 0: no such draws. 0053 % Note: when options_.ms.IndxAp=0, there is no effect from the values of options_.ms.options_.ms.options_.ms.options_.ms.indximf, IndxFore, 0054 % or options_.ms.apband. 0055 %options_.ms.apband = 1; % 1: error bands for A+; 0: no error bands for A+. 0056 %*** The following (impulse responses and forecasts) is used only if options_.ms.IndxAp=1 0057 %options_.ms.indximf = 1; % 1: generate draws of impulse responses; 0: no such draws (thus no effect 0058 % from options_.ms.imfband) 0059 %options_.ms.imfband = 1; % 1: error bands for impulse responses; 0: no error bands 0060 %options_.ms.indxfore = 0; % 1: generate draws of forecasts; 0: no such draws (thus no effect from options_.ms.foreband) 0061 %options_.ms.foreband = 0; % 1: error bands for out-of-sample forecasts; 0: no error bands 0062 % 0063 %options_.ms.indxgforhat = 1; % 1: plot unconditoinal forecasts; 0: no such plot 0064 rnum = nvar; % number of rows in the graph 0065 cnum = 1; % number of columns in the graph 0066 if rnum*cnum<nvar 0067 warning('rnum*cnum must be at least as large as nvar') 0068 disp('Hit any key to continue, or ctrl-c to abort') 0069 pause 0070 end 0071 %options_.ms.indxgimfhat = 1; % 1: plot ML impulse responses; 0: no plot 0072 %options_.ms.indxestima = 1; %1: ML estimation; 0: no estimation and data only 0073 % 0074 IndxNmlr = [1 0 0 0 0 0]; % imported by nmlzvar.m 0075 % Index for which normalization rule to choose 0076 % Only one of the elments in IndxNmlr can be non-zero 0077 % IndxNmlr(1): ML A distance rule (supposed to be the best) 0078 % IndxNmlr(2): ML Ahat distance rule (to approximate IndxNmlr(1)) 0079 % IndxNmlr(3): ML Euclidean distance rule (not invariant to scale) 0080 % IndxNmlr(4): Positive diagonal rule 0081 % IndxNmlr(5): Positive inv(A) diagonal rule (if ~IndxNmlr(5), no need for A0inu, 0082 % so let A0inu=[]) 0083 % IndxNmlr(6): Assigned postive rule (such as off-diagonal elements). Added 1/3/00 0084 0085 0086 %%%%---------------------------------------- 0087 % Hard conditions directly on variables 0088 % 0089 %options_.ms.indxgdls = 1; % 1: graph point forecast on variables; 0: disable 0090 nconstr1=nfqm; % number of the 1st set of constraints 0091 nconstr2=options_.forecast ; % number of the 2nd set of constraints 0092 nconstr=nconstr1+nconstr2; % q: 4 years -- 4*12 months. 0093 % When 0, no conditions directly on variables <<>> 0094 nconstr=0 ; %6*nconstr1; 0095 options_.ms.eq_ms = []; % location of MS equation; if [], all shocks 0096 PorR = [4*ones(nconstr1,1);2*ones(nconstr1,1);3*ones(nconstr1,1)]; % the variable conditioned. 1: Pcm; 3: FFR; 4: CPI 0097 PorR = [PorR;1*ones(nconstr1,1);5*ones(nconstr1,1);6*ones(nconstr1,1)]; 0098 %PorR = [3 5]; % the variable conditioned. 3: FFR; 5: CPI 0099 %PorR = 3; 0100 0101 %%%%---------------------------------------- 0102 % Conditions directly on future shocks 0103 % 0104 %options_.ms.cms = 0 % 1: condition on ms shocks; 0: disable this and "fidcnderr.m" gives 0105 % unconditional forecasts if nconstr = 0 as well; <<>> 0106 %options_.ms.ncms = 0; % number of the stance of policy; 0 if no tightening or loosening 0107 %options_.ms.eq_cms = 1; % location of MS shocks 0108 options_.ms.tlindx = 1*ones(1,options_.ms.ncms); % 1-by-options_.ms.ncms vector; 1: tightening; 0: loosen 0109 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 0110 % 1-by-options_.ms.ncms vector; cut-off point for MS shocks 0111 TLmean = zeros(1,options_.ms.ncms); 0112 % unconditional, i.e., 0 mean, for the final report in the paper 0113 if options_.ms.cms 0114 options_.ms.eq_ms = []; 0115 % At least at this point, it makes no sense to have DLS type of options_.ms.eq_ms; 10/12/98 0116 if all(isfinite(options_.ms.tlnumber)) 0117 for k=1:options_.ms.ncms 0118 TLmean(k) = lcnmean(options_.ms.tlnumber(k),options_.ms.tlindx(k)); 0119 % shock mean magnitude. 1: tight; 0: loose 0120 % Never used for any subsequent computation but 0121 % simply used for the final report in the paper. 0122 %options_.ms.tlnumber(k) = fzero('lcutoff',0,[],[],TLmean(k)) 0123 % get an idea about the cutoff point given TLmean instead 0124 0125 end 0126 end 0127 else 0128 options_.ms.ncms = 0; % only for the use of the graph by msprobg.m 0129 options_.ms.tlnumber = NaN*ones(1,options_.ms.ncms); 0130 % -infinity, only for the use of the graph by msprobg.m 0131 end 0132 0133 0134 %%%%---------------------------------------- 0135 % Soft conditions on variables 0136 % 0137 %cnum = 0 % # of band condtions; when 0, disable this option 0138 % Note (different from "fidencon") that each condition corres. to variable 0139 %options_.ms.banact = 1; % 1: use infor on actual; 0: preset without infor on actual 0140 if cnum 0141 banindx = cell(cnum,1); % index for each variable or conditon 0142 banstp = cell(cnum,1); % steps: annual in general 0143 banvar = zeros(cnum,1); % varables: annual in general 0144 banval = cell(cnum,1); % band value (each variable occupy a cell) 0145 badval{1} = zeros(length(banstp{1}),2); % 2: lower or higher bound 0146 0147 banstp{1} = 1:4; % 3 or 4 years 0148 banvar(1) = 3; % 3: FFR; 5: CPI 0149 if ~options_.ms.banact 0150 for i=1:length(banstp{1}) 0151 banval{1}(i,:) = [5.0 10.0]; 0152 end 0153 end 0154 end 0155 % 0156 pause(1) 0157 disp(' ') 0158 disp('For uncondtional forecasts, set nconstr = options_.ms.cms = cnum = 0') 0159 pause(1) 0160 % 0161 %================================================= 0162 %====== End of exordium ========================== 0163 %================================================= 0164 0165 0166 0167 0168 0169 %(1)-------------------------------------- 0170 % Further data analysis 0171 %(1)-------------------------------------- 0172 % 0173 if (options_.ms.freq==12) 0174 nStart=(yrStart-options_.ms.initial_year )*12+qmStart-options_.ms.initial_subperiod ; % positive number of months at the start 0175 nEnd=(yrEnd-options_.ms.final_year )*12+qmEnd-options_.ms.final_subperiod ; % negative number of months towards end 0176 elseif (options_.ms.freq==4) 0177 nStart=(yrStart-options_.ms.initial_year )*4+qmStart-options_.ms.initial_subperiod ; % positive number of months at the start 0178 nEnd=(yrEnd-options_.ms.final_year )*4+qmEnd-options_.ms.final_subperiod ; % negative number of months towards end 0179 elseif (options_.ms.freq==1) 0180 nStart=(yrStart-options_.ms.initial_year )*1+qmStart-options_.ms.initial_subperiod ; % positive number of months at the start 0181 nEnd=(yrEnd-options_.ms.final_year )*1+qmEnd-options_.ms.final_subperiod ; % negative number of months towards end 0182 else 0183 error('Error: this code is only good for monthly/quarterly/yearly data!!!') 0184 return 0185 end 0186 % 0187 if nEnd>0 || nStart<0 0188 disp('Warning: this particular sample consider is out of bounds of the data!!!') 0189 return 0190 end 0191 %*** Note, both xdgel and xdata have the same start with the specific sample 0192 xdgel=options_.data(nStart+1:nData+nEnd,options_.ms.vlist); 0193 % gel: general options_.data within sample (nSample) 0194 if ~(nSample==size(xdgel,1)) 0195 warning('The sample size (including options_.ms.nlags ) and data are incompatible') 0196 disp('Check to make sure nSample and size(xdgel,1) are the same') 0197 return 0198 end 0199 % 0200 baddata = find(isnan(xdgel)); 0201 if ~isempty(baddata) 0202 warning('Some data for this selected sample are actually unavailable.') 0203 disp('Hit any key to continue, or ctrl-c to abort') 0204 pause 0205 end 0206 % 0207 if options_.ms.initial_subperiod ==1 0208 yrB = options_.ms.initial_year ; qmB = options_.ms.initial_subperiod ; 0209 else 0210 yrB = options_.ms.initial_year +1; qmB = 1; 0211 end 0212 yrF = options_.ms.final_year ; qmF = options_.ms.final_subperiod ; 0213 [Mdate,tmp] = fn_calyrqm(options_.ms.freq,[options_.ms.initial_year options_.ms.initial_subperiod ],[options_.ms.final_year options_.ms.final_subperiod ]); 0214 xdatae=[Mdate options_.data(1:nData,options_.ms.vlist)]; 0215 % beyond sample into forecast horizon until the end of the data options_.ms.final_year :options_.ms.final_subperiod 0216 % Note: may contain NaN data. So must be careful about its use 0217 0218 %=========== Obtain prior-period, period-to-last period, and annual growth rates 0219 [yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(xdatae,options_.ms.freq,options_.ms.log_var,options_.ms.percent_var,[yrB qmB],[yrF qmF]); 0220 qdates = zeros(size(yactqmyge,1),1); 0221 for ki=1:length(qdates) 0222 qdates(ki) = yactqmyge(1,1) + (yactqmyge(1,2)+ki-2)/options_.ms.freq; 0223 end 0224 for ki=1:nvar 0225 figure 0226 plot(qdates, yactqmyge(:,2+ki)/100) 0227 xlabel(options_.ms.varlist{ki}) 0228 end 0229 save outactqmygdata.prn yactqmyge -ascii 0230 0231 0232 0233 %=========== Write the output on the screen or to a file in an organized way ============== 0234 %disp([sprintf('%4.0f %2.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',yactyrge')]) 0235 spstr1 = 'disp([sprintf('; 0236 spstr2 = '%4.0f %2.0f'; 0237 yactyrget=yactyrge'; 0238 for ki=1:length(options_.ms.vlist) 0239 if ki==length(options_.ms.vlist) 0240 spstr2 = [spstr2 ' %8.3f\n']; 0241 else 0242 spstr2 = [spstr2 ' %8.3f']; 0243 end 0244 end 0245 spstr = [spstr1 'spstr2' ', yactyrget)])']; 0246 eval(spstr) 0247 0248 % 0249 fid = fopen('outyrqm.prn','w'); 0250 %fprintf(fid,'%4.0f %2.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',yactyrge'); 0251 fpstr1 = 'fprintf(fid,'; 0252 fpstr2 = '%4.0f %2.0f'; 0253 for ki=1:nvar 0254 if ki==nvar 0255 fpstr2 = [fpstr2 ' %8.3f\n']; 0256 else 0257 fpstr2 = [fpstr2 ' %8.3f']; 0258 end 0259 end 0260 fpstr = [fpstr1 'fpstr2' ', yactyrget);']; 0261 eval(fpstr) 0262 fclose(fid); 0263 0264 0265 0266 if options_.ms.indxestima 0267 %(2)---------------------------------------------------------------------------- 0268 % Estimation 0269 % ML forecast and impulse responses 0270 % Hard or soft conditions for conditional forecasts 0271 %(2)---------------------------------------------------------------------------- 0272 % 0273 %* Arranged data information, WITHOUT dummy obs when 0 after mu is used. See fn_rnrprior_covres_dobs.m for using the dummy 0274 % observations as part of an explicit prior. 0275 [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh] = fn_dataxy(nvar,options_.ms.nlags ,xdgel,mu,0,nexo); 0276 if qmStart+options_.ms.nlags -options_.ms.dummy_obs >0 0277 qmStartEsti = rem(qmStart+options_.ms.nlags -options_.ms.dummy_obs ,options_.ms.freq); % dummy observations are included in the sample. 0278 if (~qmStartEsti) 0279 qmStartEsti = options_.ms.freq; 0280 end 0281 yrStartEsti = yrStart + floor((qmStart+options_.ms.nlags -options_.ms.dummy_obs )/(options_.ms.freq+0.01)); 0282 % + 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. 0283 else 0284 qmStartEsti = options_.ms.freq + rem(qmStart+options_.ms.nlags -options_.ms.dummy_obs ,options_.ms.freq); % dummy observations are included in the sample. 0285 if (qmStart+options_.ms.nlags -options_.ms.dummy_obs ==0) 0286 yrStartEsti = yrStart - 1; % one year back. 0287 else 0288 yrStartEsti = yrStart + floor((qmStart+options_.ms.nlags -options_.ms.dummy_obs )/(options_.ms.freq-0.01)); 0289 % - 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. 0290 end 0291 end 0292 dateswd = fn_dataext([yrStartEsti qmStartEsti],[yrEnd qmEnd],xdatae(:,[1:2])); % dates with dummies 0293 phie = [dateswd phi]; 0294 ye = [dateswd y]; 0295 0296 %* Obtain linear restrictions 0297 [Uiconst,Viconst,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,nvar,nexo,options_.ms ); 0298 if min(n0)==0 0299 disp(' ') 0300 warning('A0: restrictions in dlrprior.m give no free parameter in one of equations') 0301 disp('Press ctrl-c to abort') 0302 pause 0303 elseif min(np)==0 0304 disp(' ') 0305 warning('Ap: Restrictions in dlrprior.m give no free parameter in one of equations') 0306 disp('Press ctrl-c to abort') 0307 pause 0308 end 0309 0310 if options_.ms.contemp_reduced_form 0311 Uiconst=cell(nvar,1); Viconst=cell(ncoef,1); 0312 for kj=1:nvar 0313 Uiconst{kj} = eye(nvar); Viconst{kj} = eye(ncoef); 0314 end 0315 end 0316 0317 if options_.ms.bayesian_prior 0318 %*** Obtains asymmetric prior (with no linear restrictions) with dummy observations as part of an explicit prior (i.e, 0319 % reflected in Hpmulti and Hpinvmulti). See Forecast II, pp.69a-69b for details. 0320 if 1 % Liquidity effect prior on both MS and MD equations. 0321 [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior_covres_dobs(nvar,options_.ms.freq,options_.ms.nlags ,xdgel,mu,indxDummy,hpmsmd,indxmsmdeqn); 0322 else 0323 [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] = fn_rnrprior(nvar,options_.ms.freq,options_.ms.nlags ,xdgel,mu); 0324 end 0325 0326 %*** Combines asymmetric prior with linear restrictions 0327 [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Uiconst,Viconst,Pi,H0multi,Hpmulti,nvar); 0328 0329 %*** Obtains the posterior matrices for estimation and inference 0330 [Pmat,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0invtld,Hpinvtld,Uiconst,Viconst); 0331 0332 if options_.ms.contemp_reduced_form 0333 %*** Obtain the ML estimate 0334 A0hatinv = chol(H0inv{1}/fss); % upper triangular but lower triangular choleski 0335 A0hat=inv(A0hatinv); 0336 a0indx = find(A0hat); 0337 else 0338 %*** Obtain the ML estimate 0339 % load idenml 0340 x = 10*rand(sum(n0),1); 0341 H0 = eye(sum(n0)); 0342 crit = 1.0e-9; 0343 nit = 10000; 0344 % 0345 [fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = ... 0346 csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Uiconst,nvar,n0,fss,H0inv); 0347 0348 A0hat = fn_tran_b2a(xhat,Uiconst,nvar,n0) 0349 A0hatinv = inv(A0hat); 0350 fhat 0351 xhat 0352 grad 0353 itct 0354 fcount 0355 retcodehat 0356 save outm.mat xhat A0hat A0hatinv grad fhat itct itct fcount retcodehat 0357 end 0358 else 0359 %*** Obtain the posterior matrices for estimation and inference 0360 [Pmat,H0inv,Hpinv] = fn_dlrpostr(xtx,xty,yty,Uiconst,Viconst); 0361 0362 if options_.ms.contemp_reduced_form 0363 %*** Obtain the ML estimate 0364 A0hatinv = chol(H0inv{1}/fss); % upper triangular but lower triangular choleski 0365 A0hat=inv(A0hatinv); 0366 a0indx = find(A0hat); 0367 else 0368 %*** Obtain the ML estimate 0369 % load idenml 0370 x = 10*rand(sum(n0),1); 0371 H0 = eye(sum(n0)); 0372 crit = 1.0e-9; 0373 nit = 10000; 0374 % 0375 [fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = ... 0376 csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Uiconst,nvar,n0,fss,H0inv); 0377 0378 A0hat = fn_tran_b2a(xhat,Uiconst,nvar,n0) 0379 A0hatinv = inv(A0hat); 0380 fhat 0381 xhat 0382 grad 0383 itct 0384 fcount 0385 retcodehat 0386 save outm.mat xhat A0hat A0hatinv grad fhat itct itct fcount retcodehat 0387 end 0388 end 0389 0390 %**** impulse responses 0391 swish = A0hatinv; % each column corresponds to an equation 0392 if options_.ms.contemp_reduced_form 0393 xhat = A0hat(a0indx); 0394 Bhat=Pmat{1}; 0395 Fhat = Bhat*A0hat 0396 ghat = NaN; 0397 else 0398 xhat = fn_tran_a2b(A0hat,Uiconst,nvar,n0); 0399 [Fhat,ghat] = fn_gfmean(xhat,Pmat,Viconst,nvar,ncoef,n0,np); 0400 if options_.ms.cross_restrictions 0401 Fhatur0P = Fhat; % ur: unrestriced across A0 and A+ 0402 for ki = 1:size(ixmC0Pres,1) % loop through the number of equations in which 0403 % cross-A0-A+ restrictions occur. See St. Louis Note p.5. 0404 ixeq = ixmC0Pres{ki}(1,1); % index for the jth equation in consideration. 0405 Lit = Viconst{ixeq}(ixmC0Pres{ki}(:,2),:); % transposed restriction matrix Li 0406 % V_j(i,:) in f_j(i) = V_j(i,:)*g_j 0407 ci = ixmC0Pres{ki}(:,4) .* A0hat(ixmC0Pres{ki}(:,3),ixeq); 0408 % s * a_j(h) in the restriction f_j(i) = s * a_j(h). 0409 LtH = Lit/Hpinv{ixeq}; 0410 HLV = LtH'/(LtH*Lit'); 0411 gihat = Viconst{ixeq}'*Fhatur0P(:,ixeq); 0412 Fhat(:,ixeq) = Viconst{ixeq}*(gihat + HLV*(ci-Lit*gihat)); 0413 end 0414 end 0415 Fhat 0416 Bhat = Fhat/A0hat; % ncoef-by-nvar reduced form lagged parameters. 0417 end 0418 nn = [nvar options_.ms.nlags imstp]; 0419 imfhat = fn_impulse(Bhat,swish,nn); % in the form that is congenial to RATS 0420 imf3hat=reshape(imfhat,size(imfhat,1),nvar,nvar); 0421 % imf3: row--steps, column--nvar responses, 3rd dimension--nvar shocks 0422 imf3shat=permute(imf3hat,[1 3 2]); 0423 % imf3s: permuted so that row--steps, column--nvar shocks, 0424 % 3rd dimension--nvar responses 0425 % Note: reshape(imf3s(1,:,:),nvar,nvar) = A0in (columns -- equations) 0426 if options_.ms.indxgimfhat 0427 figure 0428 end 0429 scaleout = fn_imcgraph(imfhat,nvar,imstp,xlab,ylab,options_.ms.indxgimfhat); 0430 imfstd = max(abs(scaleout)'); % row: nvar (largest number); used for standard deviations 0431 0432 % 0433 % %**** save stds. of both data and impulse responses in idfile1 0434 % temp = [std(yactqmyge(:,3:end)); std(yactyrge(:,3:end)); imfstd]; %<<>> 0435 % save idenyimstd.prn temp -ascii % export forecast and impulse response to the file "idenyimstd.prn", 3-by-nvar 0436 % % 0437 % %**** save stds. of both data and impulse responses in idfile1 0438 % temp = [std(yactqmyge(:,3:end)); std(yactyrge(:,3:end)); imfstd]; %<<>> 0439 % save idenyimstd.prn temp -ascii % export forecast and impulse response to the file "idenyimstd.prn", 3-by-nvar 0440 % if options_.ms.indxparr 0441 % idfile1='idenyimstd'; 0442 % end 0443 0444 %===================================== 0445 % Now, out-of-sample forecasts. Note: Hm1t does not change with A0. 0446 %===================================== 0447 % 0448 % * updating the last row of X (phi) with the current (last row of) y. 0449 tcwx = nvar*options_.ms.nlags ; % total coefficeint without exogenous variables 0450 phil = phi(size(phi,1),:); 0451 phil(nvar+1:tcwx) = phil(1:tcwx-nvar); 0452 phil(1:nvar) = y(end,:); 0453 %*** exogenous variables excluding constant terms 0454 if (nexo>1) 0455 Xexoe = fn_dataext([yrEnd qmEnd],[yrEnd qmEnd],xdatae(:,[1:2 2+nvar+1:2+nvar+nexo-1])); 0456 phil(1,tcwx+1:tcwx+nexo-1) = Xexoe(1,3:end); 0457 end 0458 % 0459 %*** ML unconditional point forecast 0460 nn = [nvar options_.ms.nlags nfqm]; 0461 if nexo<2 0462 yforehat = fn_forecast(Bhat,phil,nn); % nfqm-by-nvar, in log 0463 else 0464 Xfexoe = fn_dataext(fdates(1,:),fdates(numel(fdates),:),xdatae(:,[1:2 2+nvar+1:2+nvar+nexo-1])); 0465 %Xfexoe = fn_dataext(fdates(1,:),fdates(end,:),xdatae(:,[1:2 2+nvar+1:2+nvar+nexo-1])); 0466 yforehat = fn_forecast(Bhat,phil,nn,nexo,Xfexoe(:,3:end)); % nfqm-by-nvar, in log 0467 end 0468 yforehate = [fdates yforehat]; 0469 % 0470 yact1e = fn_dataext([yrEnd-nayr 1],[yrEnd qmEnd],xdatae(:,1:nvar+2)); 0471 if options_.ms.real_pseudo_forecast 0472 %yact2e = fn_dataext([yrEnd-nayr 1],E2yrqm,xdatae); 0473 yact2e = fn_dataext([yrEnd-nayr 1],[fdates(end,1) options_.ms.freq],xdatae(:,1:nvar+2)); 0474 else 0475 yact2e=yact1e; 0476 end 0477 yafhate = [yact1e; yforehate]; % actual and forecast 0478 % 0479 %===== Converted to mg, qg, and calendar yg 0480 % 0481 [yafyrghate,yafyrhate,yafqmyghate] = fn_datana(yafhate,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno)); 0482 % actual and forecast growth rates 0483 [yact2yrge,yact2yre,yact2qmyge] = fn_datana(yact2e,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno)); 0484 % only actual growth rates 0485 yafyrghate 0486 if options_.ms.indxgforhat 0487 keyindx = [1:nvar]; 0488 conlab=['unconditional']; 0489 0490 figure 0491 yafyrghate(:,3:end) = yafyrghate(:,3:end)/100; 0492 yact2yrge(:,3:end) = yact2yrge(:,3:end)/100; 0493 fn_foregraph(yafyrghate,yact2yrge,keyindx,rnum,cnum,options_.ms.freq,ylab,forelabel,conlab) 0494 end 0495 0496 %------------------------------------------------- 0497 % Setup for point conditional forecast 0498 % ML Conditional Forecast 0499 %------------------------------------------------- 0500 % 0501 %% See Zha's note "Forecast (1)" p. 5, RATS manual (some errors in RATS), etc. 0502 % 0503 %% Some notations: y(t+1) = y(t)B1 + e(t+1)inv(A0). e(t+1) is 1-by-n. 0504 %% Let r(t+1)=e(t+1)inv(A0) + e(t+2)C + .... where inv(A0) is impulse 0505 %% response at t=1, C at t=2, etc. The row of inv(A0) or C is 0506 %% all responses to one shock. 0507 %% Let r be q-by-1 (such as r(1) = r(t+1) 0508 %% = y(t+1) (constrained) - y(t+1) (forecast)). 0509 %% Use impulse responses to find out R (k-by-q) where k=nvar*nsteps 0510 %% where nsteps the largest constrained step. The key of the program 0511 %% is to creat R using impulse responses 0512 %% Optimal solution for shock e where R'*e=r and e is k-by-1 is 0513 %% e = R*inv(R'*R)*r. 0514 % 0515 0516 if (nconstr > 0) 0517 %*** initializing 0518 stepcon=cell(nconstr,1); % initializing, value y conditioned 0519 valuecon=zeros(nconstr,1); % initializing, value y conditioned 0520 varcon=zeros(nconstr,1); % initializing, endogous variables conditioned 0521 varcon(:)=PorR; % 1: Pcm; 3: FFR; 5: CPI 0522 0523 % 0524 for i=1:nconstr 0525 if i<=nconstr1 0526 stepcon{i}=i; % FFR 0527 elseif i<=2*nconstr1 0528 stepcon{i}=i-nconstr1; % FFR 0529 elseif i<=3*nconstr1 0530 stepcon{i}=i-2*nconstr1; % FFR 0531 elseif i<=4*nconstr1 0532 stepcon{i}=i-3*nconstr1; % FFR 0533 elseif i<=5*nconstr1 0534 stepcon{i}=i-4*nconstr1; % FFR 0535 elseif i<=6*nconstr1 0536 stepcon{i}=i-5*nconstr1; % FFR 0537 end 0538 end 0539 0540 % for i=1:nconstr 0541 % stepcon{i}=i; % FFR 0542 % end 0543 0544 % bend=12; 0545 % stepcon{1}=[1:bend]'; % average over 0546 % stepcon{nconstr1+1}=[1:options_.ms.freq-qmSub]'; % average over the remaing months in 1st forecast year 0547 % stepcon{nconstr1+2}=[options_.ms.freq-qmSub+1:options_.ms.freq-qmSub+12]'; % average over 12 months next year 0548 % stepcon{nconstr1+3}=[options_.ms.freq-qmSub+13:options_.ms.freq-qmSub+24]'; % average over 12 months. 3rd year 0549 % stepcon{nconstr1+4}=[options_.ms.freq-qmSub+25:options_.ms.freq-qmSub+36]'; % average over 12 months. 4th year 0550 0551 % %**** avearage condition over, say, options_.ms.freq periods 0552 % if qmEnd==options_.ms.freq 0553 % stepcon{1}=[1:options_.ms.freq]'; % average over the remaing periods in 1st forecast year 0554 % else 0555 % stepcon{1}=[1:options_.ms.freq-qmEnd]'; % average over the remaing periods in 1st forecast year 0556 % end 0557 % for kj=2:nconstr 0558 % stepcon{kj}=[length(stepcon{kj-1})+1:length(stepcon{kj-1})+options_.ms.freq]'; % average over 12 months next year 0559 % end 0560 0561 if options_.ms.real_pseudo_forecast 0562 % %*** conditions in every period 0563 % for i=1:nconstr 0564 % valuecon(i) = yact(actup+i,varcon(i)); 0565 % %valuecon(i) = mean( yact(actup+1:actup+bend,varcon(i)) ); 0566 % %valuecon(i) = 0.060; % 95:01 0567 % %valuecon(i) = (0.0475+0.055)/2; % 94:10 0568 % end 0569 0570 % %*** average condtions over,say, options_.ms.freq periods. 0571 % for i=nconstr1+1:nconstr1+nconstr2 0572 % i=1; 0573 % valuecon(nconstr1+i) = ( ( mean(ylast12Cal(:,varcon(nconstr1+i)),1) + ... 0574 % log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100) )*options_.ms.freq - ... 0575 % yCal_1(:,varcon(nconstr1+i)) ) ./ length(stepcon{nconstr1+i}); 0576 % % the same as unconditional "yactCalyg" 1st calendar year 0577 % i=2; 0578 % valuecon(nconstr1+i) = mean(ylast12Cal(:,varcon(nconstr1+i))) + ... 0579 % log(1+yactCalyg(yAg-yFg+1,varcon(nconstr1+i))/100) ... 0580 % + log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100); 0581 % % the same as actual "yactCalgy" 2nd calendar year 0582 % i=3; 0583 % valuecon(nconstr1+i) = valuecon(nconstr1+i-1) + ... 0584 % log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100); 0585 % % the same as actual "yactCalgy" 3rd calendar year 0586 % %i=4; 0587 % %valuecon(nconstr1+i) = valuecon(nconstr1+i-1) + ... 0588 % % log(1+yactCalyg(yAg-yFg+i,varcon(nconstr1+i))/100); 0589 % % the same as actual "yactCalgy" 4th calendar year 0590 % end 0591 0592 %*** conditions in every period 0593 vpntM = fn_dataext(E1yrqm, E2yrqm,xdatae); % point value matrix with dates 0594 % vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates 0595 for i=1:nconstr 0596 if i<=nconstr1 0597 valuecon(i) = vpntM(i,2+varcon(i)); % 2: first 2 elements are dates 0598 elseif i<=2*nconstr1 0599 valuecon(i) = vpntM(i-nconstr1,2+varcon(i)); 0600 elseif i<=3*nconstr1 0601 valuecon(i) = vpntM(i-2*nconstr1,2+varcon(i)); 0602 elseif i<=4*nconstr1 0603 valuecon(i) = vpntM(i-3*nconstr1,2+varcon(i)); 0604 elseif i<=5*nconstr1 0605 valuecon(i) = vpntM(i-4*nconstr1,2+varcon(i)); 0606 elseif i<=6*nconstr1 0607 valuecon(i) = vpntM(i-5*nconstr1,2+varcon(i)); 0608 end 0609 end 0610 0611 % %*** average condtions over,say, options_.ms.freq periods. 0612 % if qmEnd==options_.ms.freq 0613 % vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates 0614 % valuecon(1) = vaveM(1,2+varcon(1)); % 2: first 2 elements are dates 0615 % else 0616 % vaveM = fn_dataext([yrEnd 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates 0617 % yactrem = fn_dataext([yrEnd qmEnd+1],[yrEnd options_.ms.freq],xdatae); 0618 % valuecon(1) = sum(yactrem(:,2+varcon(1)),1)/length(stepcon{1}); 0619 % % 2: first 2 elements are dates 0620 % end 0621 % for kj=2:nconstr 0622 % valuecon(kj) = vaveM(kj,2+varcon(kj)); % 2: first 2 elements are dates 0623 % end 0624 else 0625 vpntM = dataext([yrEnd qmEnd+1],[yrEnd qmEnd+2],xdatae); % point value matrix with dates 0626 for i=1:nconstr 0627 if i<=nconstr1 0628 valuecon(i) = vpntM(i,2+varcon(i)); % 2: first 2 elements are dates; Poil 0629 elseif i<=2*nconstr1 0630 valuecon(i) = vpntM(i-nconstr1,2+varcon(i)); % 2: first 2 elements are dates; M2 0631 elseif i<=3*nconstr1 0632 valuecon(i) = vpntM(i-2*nconstr1,2+varcon(i)); % 2: first 2 elements are dates; FFR 0633 elseif i<=4*nconstr1 0634 valuecon(i) = vpntM(i-3*nconstr1,2+varcon(i)); % 2: first 2 elements are dates; CPI 0635 elseif i<=5*nconstr1 0636 valuecon(i) = vpntM(i-4*nconstr1,2+varcon(i)); % 2: first 2 elements are dates; U 0637 elseif i<=5*nconstr1+nconstr2 0638 valuecon(i)=xdata(end,5)+(i-5*nconstr1)*log(1.001)/options_.ms.freq; %CPI 0639 elseif i<=5*nconstr1+2*nconstr2 0640 valuecon(i)=0.0725; %FFR 0641 else 0642 valuecon(i)=xdata(end,6)+(i-5*nconstr1-2*nconstr2)*0.01/nfqm; %U 0643 end 0644 end 0645 %valuecon(i) = 0.060; % 95:01 0646 end 0647 else 0648 valuecon = []; 0649 stepcon = []; 0650 varcon = []; 0651 end 0652 0653 nstepsm = 0; % initializing, the maximum step in all constraints 0654 for i=1:nconstr 0655 nstepsm = max([nstepsm max(stepcon{i})]); 0656 end 0657 0658 if cnum 0659 if options_.ms.real_pseudo_forecast && options_.ms.banact 0660 for i=1:length(banstp{1}) 0661 banval{1}(1:length(banstp{1}),1) = ... 0662 yactCalyg(yAg-yFg+1:yAg-yFg+length(banstp{1}),banvar(1)) - 2; 0663 banval{1}(1:length(banstp{1}),2) = ... 0664 yactCalyg(yAg-yFg+1:yAg-yFg+length(banstp{1}),banvar(1)) + 2; 0665 end 0666 end 0667 end 0668 0669 0670 %=================================================== 0671 % ML conditional forecast 0672 %=================================================== 0673 %/* 0674 [ychat,Estr,rcon] = fn_fcstidcnd(valuecon,stepcon,varcon,nstepsm,... 0675 nconstr,options_.ms.eq_ms,nvar,options_.ms.nlags ,phil,0,0,yforehat,imf3shat,A0hat,Bhat,... 0676 nfqm,options_.ms.tlindx,options_.ms.tlnumber,options_.ms.ncms,options_.ms.eq_cms); 0677 ychate = [fdates ychat]; 0678 yachate = [yact1e; ychate]; % actual and condtional forecast 0679 %===== Converted to mg, qg, and calendar yg 0680 [yacyrghate,yacyrhate,yacqmyghate] = fn_datana(yachate,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno)); 0681 % actual and conditional forecast growth rates 0682 if options_.ms.indxgdls && nconstr 0683 keyindx = [1:nvar]; 0684 % conlab=['conditional on' ylab{PorR(1)}]; 0685 conlab=['v-conditions']; 0686 0687 figure 0688 fn_foregraph(yafyrghate,yact2yrge,keyindx,rnum,cnum,options_.ms.freq,ylab,forelabel,conlab) 0689 end 0690 0691 if options_.ms.ncsk 0692 Estr = zeros(nfqm,nvar); 0693 Estr(1:2,:) = [ 0694 -2.1838 -1.5779 0.53064 -0.099425 -0.69269 -1.0391 0695 1.9407 3.3138 -0.10563 -0.55457 -0.68772 1.3534 0696 ]; 0697 Estr(3:6,3) = [0.5*ones(1,4)]'; % MD shocks 0698 0699 Estr(3:10,2) = [1.5 1.5 1.5*ones(1,6)]'; % MS shocks 0700 0701 %Estr(3:6,6) = 1*ones(4,1); % U shocks 0702 %Estr(8:11,4) = 1*ones(4,1); % y shocks 0703 0704 %Estr(3:10,2) = [2.5 2.5 1.5*ones(1,6)]'; % MS shocks alone 0705 0706 nn = [nvar options_.ms.noptions_.ms.nlags nfqm]; 0707 ycEhat = forefixe(A0hat,Bhat,phil,nn,Estr); 0708 ycEhate = [fdates ycEhat]; 0709 yacEhate = [yact1e; ycEhate]; % actual and condtional forecast 0710 %===== Converted to mg, qg, and calendar yg 0711 [yacEyrghate,yacEyrhate,yacEqmyghate] = datana(yacEhate,options_.ms.freq,options_.ms.log_var(1:nlogeno),options_.ms.percent_var(1:npereno)); 0712 % actual and conditional forecast growth rates 0713 disp([sprintf('%4.0f %2.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',yacEyrghate')]) 0714 0715 if 1 0716 keyindx = [1:nvar]; 0717 % conlab=['conditional on' ylab{PorR(1)}]; 0718 conlab=['shock-conditions']; 0719 0720 figure 0721 gyrfore(yacEyrghate,yact2yrge,keyindx,rnum,cnum,ylab,forelabel,conlab) 0722 end 0723 end 0724 0725 %----------------------------------------------------------- 0726 % Compute structural shocks for the whole sample period excluding dummy observations. 0727 %----------------------------------------------------------- 0728 ywod = y(options_.ms.dummy_obs +1:end,:); % without dummy observations 0729 phiwod=phi(options_.ms.dummy_obs +1:end,:); % without dummy observations 0730 eplhat=ywod*A0hat-phiwod*Fhat; 0731 qmStartWod = mod(qmStart+options_.ms.nlags ,options_.ms.freq); 0732 if (~qmStartWod) 0733 qmStartWod = options_.ms.freq; 0734 end 0735 yrStartWod = yrStart + floor((qmStart+options_.ms.nlags -1)/options_.ms.freq); 0736 dateswod = fn_dataext([yrStartWod qmStartWod],[yrEnd qmEnd],xdatae(:,[1:2])); 0737 eplhate = [dateswod eplhat]; 0738 0739 Aphat = Fhat; 0740 end 0741 0742 %---------------------------------------- 0743 % Tests for LR, HQ, Akaike, Schwarz. The following gives a guidance. 0744 % But the computation has to be done in a different M file by exporting fhat's 0745 % from different idfile's. 0746 %---------------------------------------- 0747 % 0748 %if ~options_.ms.contemp_reduced_form 0749 % SpHR=A0in'*A0in; 0750 %end 0751 %% 0752 %if ~isnan(SpHR) && ~options_.ms.contemp_reduced_form 0753 % warning(' ') 0754 % disp('Make sure you run the program with options_.ms.contemp_reduced_form =1 first.') 0755 % disp('Otherwise, the following test results such as Schwartz are incorrect.') 0756 % disp('All other results such as A0ml and imfs, however, are correct.') 0757 % disp('Press anykey to contintue or ctrl-c to stop now') 0758 % pause 0759 0760 % load SpHUout 0761 0762 % logLHU=-fss*sum(log(diag(chol(SpHU)))) -0.5*fss*nvar % unrestricted logLH 0763 0764 % logLHR=-fhat % restricted logLH 0765 % tra = reshape(SpHU,nvar*nvar,1)'*reshape(A0*A0',nvar*nvar,1); 0766 % df=(nvar*(nvar+1)/2 - length(a0indx)); 0767 % S=2*(logLHU-logLHR); 0768 % SC = (nvar*(nvar+1)/2 - length(a0indx)) * log(fss); 0769 % disp(['T -- effective sample size: ' num2str(fss)]) 0770 % disp(['Trace in the overidentified posterior: ' num2str(tra)]) 0771 % disp(['Chi2 term -- 2*(logLHU-logLHR): ' num2str(S)]) 0772 % disp(['Degrees of freedom: ' num2str(df)]) 0773 % disp(['SC -- df*log(T): ' num2str(SC)]) 0774 % disp(['Akaike -- 2*df: ' num2str(2*df)]) 0775 % disp(['Classical Asymptotic Prob at chi2 term: ' num2str(cdf('chi2',S,df))]) 0776 0777 % %*** The following is the eigenanalysis in the difference between 0778 % %*** unrestricted (U) and restricted (R) 0779 % norm(A0'*SpHU*A0-diag(diag(ones(6))))/6; 0780 % norm(SpHU-A0in'*A0in)/6; 0781 0782 % corU = corr(SpHU); 0783 % corR = corr(SpHR); 0784 0785 % [vU,dU]=eigsort(SpHU,1); 0786 % [vR,dR]=eigsort(SpHR,1); 0787 0788 % [log(diag(dU)) log(diag(dR)) log(diag(dU))-log(diag(dR))]; 0789 0790 % sum(log(diag(dU))); 0791 % sum(log(diag(dR))); 0792 %else 0793 % disp('To run SC test, turn options_.ms.contemp_reduced_form =1 first and then turn options_.ms.contemp_reduced_form =0') 0794 %end 0795 0796 0797 %***** Simply regression 0798 %X=[phi(:,3) y(:,2)-phi(:,2) y(:,1)-phi(:,7) ones(fss,1)]; 0799 %� Y=y(:,3); 0800 %� b=regress(Y,X) 0801 0802 %=== Computes the roots for the whole system. 0803 rootsinv = fn_varoots(Bhat,nvar,options_.ms.nlags ) 0804 abs(rootsinv) 0805 0806 0807 bhat =xhat; 0808 n0const=n0; % For constant parameter models. 0809 n0const=n0; % For constant parameter models. 0810 npconst=np; % For constant parameter models. 0811 save outdata_a0dp_const.mat A0hat bhat Aphat n0const