1function [Pi_bar,H0tldcell_inv,Hptldcell_inv] ... 2 = fn_rnrprior_covres_dobs_tv(nvar,nStates,q_m,lags,xdgel,mu,indxDummy,Ui,Vi,hpmsmd,indxmsmdeqn,nexo,asym0,asymp) 3% Differs from fn_rnrprior_covres_dobs_tv(): no linear restrictions (Ui and Vi) have applied yet to this function, but 4% linear restrictions are incorported in fn_rnrprior_covres_dobs_tv(). 5% 6% Only works for the nexo=1 (constant term) case. To extend this to other exogenous variables, see fn_dataxy.m. 01/14/03. 7% Differs from fn_rnrprior_covres_tv.m in that dummy observations are included as part of the explicit prior. See Forcast II, pp.68-69b. 8% Exports random Bayesian prior of Sims and Zha with asymmetric rior with linear restrictions already applied 9% and with dummy observations (i.e., mu(5) and mu(6)) used as part of an explicit prior. 10% This function allows for prior covariances for the MS and MD equations to achieve liquidity effects. 11% See Waggoner and Zha's Gibbs sampling paper and TVBVAR NOTES pp. 71k.0 and 50-61. 12% 13% nvar: number of endogenous variables 14% nStates: Number of states. 15% q_m: quarter or month 16% lags: the maximum length of lag 17% xdgel: T*nvar endogenous-variable matrix of raw or original data (no manipulation involved) with sample size including lags. 18% Order of columns: (1) nvar endogenous variables; (2) constants will be automatically put in the last column. 19% Used only to get variances of residuals for mu(1)-mu(5) and for dummy observations mu(5) and mu(6). 20% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where 21% mu(5) and mu(6) are NOT used here. See fn_dataxy.m for using mu(5) and mu(6). 22% mu(1): overall tightness and also for A0; (0.57) 23% mu(2): relative tightness for A+; (0.13) 24% mu(3): relative tightness for the constant term; (0.1). NOTE: for other 25% exogenous terms, the variance of each exogenous term must be taken into 26% acount to eliminate the scaling factor. 27% mu(4): tightness on lag decay; (1) 28% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5) 29% mu(6): weight on single dummy initial observation including constant 30% (cointegration, unit roots, and stationarity); (5) 31% NOTE: for this function, mu(5) and mu(6) are not used. See fn_dataxy.m for using mu(5) and mu(6). 32% indxDummy: 1: uses dummy observations to form part of an explicit prior; 0: no dummy observations as part of the prior. 33% Ui: nvar-by-1 cell. In each cell, nvar-by-qi*si orthonormal basis for the null of the ith 34% equation contemporaneous restriction matrix where qi is the number of free parameters 35% within the state and si is the number of free states. 36% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector 37% of total original parameters and bi is a vector of free parameters. When no 38% restrictions are imposed, we have Ui = I. There must be at least one free 39% parameter left for the ith equation in the order of [a_i for 1st state, ..., a_i for last state]. 40% Vi: nvar-by-1 cell. In each cell, k-by-ri*ti orthonormal basis for the null of the ith 41% equation lagged restriction matrix where k is a total of exogenous variables and 42% ri is the number of free parameters within the state and ti is the number of free states. 43% With this transformation, we have fi = Vi*gi 44% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a 45% vector of free parameters. The ith equation is in the order of [nvar variables 46% for 1st lag and 1st state, ..., nvar variables for last lag and 1st state, const for 1st state, nvar 47% variables for 1st lag and 2nd state, nvar variables for last lag and 2nd state, const for 2nd state, and so on]. 48% hpmsmd: 2-by-1 hyperparameters with -1<h1=hpmsmd(1)<=0 for the MS equation and 0<=h2=hpmsmd(2)<1 the MD equation. Consider a1*R + a2*M. 49% The term h1*var(a1)*var(a2) is the prior covariance of a1 and a2 for MS, equivalent to penalizing the same sign of a1 and a2. 50% The term h2*var(a1)*var(a2) is the prior covariance of a1 and a2 for MD, equivalent to penalizing opposite signs of a1 and a2. 51% This will give us a liquidity effect. If hpmsmd=0, no such restrictions will be imposed. 52% indxmsmdeqn: 4-by-1 index for the locations of the MS and MD equation and for the locations of M and R. 53% indxmsmdeqn(1) for MS and indxmsmdeqn(2) for MD. 54% indxmsmdeqn(3) for M and indxmsmdeqn(4) for R. 55% nexo: number of exogenous variables (if not specified, nexo=1 (constant) by default). 56% The constant term is always put to the last of all endogenous and exogenous variables. 57% asym0: nvar-by-nvar asymmetric prior on A0. Column -- equation. 58% If ones(nvar,nvar), symmetric prior; if not, relative (asymmetric) tightness on A0. 59% asymp: ncoef-1-by-nvar asymmetric prior on A+ bar constant. Column -- equation. 60% If ones(ncoef-1,nvar), symmetric prior; if not, relative (asymmetric) tightness on A+. 61% -------------------- 62% Pi_bar: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations 63% H0tldcell_inv: cell(nvar,1). The ith cell represents the ith equation, where the dim is 64% qi*si-by-qi*si. The inverse of H0tld on p.60. 65% Hptldcell_inv: cell(nvar,1). The ith cell represents the ith equation, where the dim is 66% ri*ti-by-ri*ti.The inverse of Hptld on p.60. 67% 68% Tao Zha, February 2000. Revised, September 2000, 2001, February, May 2003. 69% 70% Copyright (C) 1997-2012 Tao Zha 71% 72% This free software: you can redistribute it and/or modify 73% it under the terms of the GNU General Public License as published by 74% the Free Software Foundation, either version 3 of the License, or 75% (at your option) any later version. 76% 77% It is distributed in the hope that it will be useful, 78% but WITHOUT ANY WARRANTY; without even the implied warranty of 79% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 80% GNU General Public License for more details. 81% 82% If you did not received a copy of the GNU General Public License 83% with this software, see <http://www.gnu.org/licenses/>. 84% 85 86 87if nargin<=11, nexo=1; end % <<>>1 88ncoef = nvar*lags+nexo; % Number of coefficients in *each* equation for each state, RHS coefficients only. 89ncoefsts = nStates*ncoef; % Number of coefficients in *each* equation in all states, RHS coefficients only. 90 91H0tldcell_inv=cell(nvar,1); % inv(H0tilde) for different equations under asymmetric prior. 92Hptldcell_inv=cell(nvar,1); % inv(H+tilde) for different equations under asymmetric prior. 93 94%*** Constructing Pi_bar for the ith equation under the random walk assumption 95Pi_bar = zeros(ncoef,nvar); % same for all equations 96Pi_bar(1:nvar,1:nvar) = eye(nvar); % random walk 97 98% 99%@@@ Prepared for Bayesian prior 100% 101% 102% ** monthly lag decay in order to match quarterly decay: a*exp(bl) where 103% ** l is the monthly lag. Suppose quarterly decay is 1/x where x=1,2,3,4. 104% ** Let the decay of l1 (a*exp(b*l1)) match that of x1 (say, beginning: 1/1) 105% ** and the decay of l2 (a*exp(b*l2)) match that of x2 (say, end: 1/5), 106% ** we can solve for a and b which are 107% ** b = (log_x1-log_x2)/(l1-l2), and a = x1*exp(-b*l1). 108if q_m==12 109 l1 = 1; % 1st month == 1st quarter 110 xx1 = 1; % 1st quarter 111 l2 = lags; % last month 112 xx2 = 1/((ceil(lags/3))^mu(4)); % last quarter 113 %xx2 = 1/6; % last quarter 114 % 3rd quarter: i.e., we intend to let decay of the 6th month match 115 % that of the 3rd quarter, so that the 6th month decays a little 116 % faster than the second quarter which is 1/2. 117 if lags==1 118 b = 0; 119 else 120 b = (log(xx1)-log(xx2))/(l1-l2); 121 end 122 a = xx1*exp(-b*l1); 123end 124 125 126 127% 128% *** specify the prior for each equation separately, SZ method, 129% ** get the residuals from univariate regressions. 130% 131sgh = zeros(nvar,1); % square root 132sgsh = sgh; % square 133nSample=size(xdgel,1); % sample size-lags 134yu = xdgel; 135C = ones(nSample,1); 136for k=1:nvar 137 [Bk,ek,junk1,junk2,junk3,junk4] = sye([yu(:,k) C],lags); 138 clear Bk junk1 junk2 junk3 junk4; 139 sgsh(k) = ek'*ek/(nSample-lags); 140 sgh(k) = sqrt(sgsh(k)); 141end 142% ** prior variance for A0(:,1), same for all equations!!! 143sg0bid = zeros(nvar,1); % Sigma0_bar diagonal only for the ith equation 144for j=1:nvar 145 sg0bid(j) = 1/sgsh(j); % sgsh = sigmai^2 146end 147% ** prior variance for lagged and exogeous variables, same for all equations 148sgpbid = zeros(ncoef,1); % Sigma_plus_bar, diagonal, for the ith equation 149for i = 1:lags 150 if (q_m==12) 151 lagdecay = a*exp(b*i*mu(4)); 152 end 153 % 154 for j = 1:nvar 155 if (q_m==12) 156 % exponential decay to match quarterly decay 157 sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation 158 elseif (q_m==4) 159 sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation 160 else 161 error('Incompatibility with lags, check the possible errors!!!') 162 %warning('Incompatibility with lags, check the possible errors!!!') 163 %return 164 end 165 end 166end 167% 168if indxDummy % Dummy observations as part of the explicit prior. 169 ndobs=nvar+1; % Number of dummy observations: nvar unit roots and 1 cointegration prior. 170 phibar = zeros(ndobs,ncoef); 171 %* constant term 172 const = ones(nvar+1,1); 173 const(1:nvar) = 0.0; 174 phibar(:,ncoef) = const; % the first nvar periods: no or zero constant! 175 176 xdgelint = mean(xdgel(1:lags,:),1); % mean of the first lags initial conditions 177 %* Dummies 178 for k=1:nvar 179 for m=1:lags 180 phibar(ndobs,nvar*(m-1)+k) = xdgelint(k); 181 phibar(k,nvar*(m-1)+k) = xdgelint(k); 182 % <<>> multiply hyperparameter later 183 end 184 end 185 phibar(1:nvar,:) = 1*mu(5)*phibar(1:nvar,:); % standard Sims and Zha prior 186 phibar(ndobs,:) = mu(6)*phibar(ndobs,:); 187 [phiq,phir]=qr(phibar,0); 188 xtxbar=phir'*phir; % phibar'*phibar. ncoef-by-ncoef. Reduced (not full) rank. See Forcast II, pp.69-69b. 189end 190 191 192 193 194 195%================================================= 196% Computing the (prior) covariance matrix for the posterior of A0, no data yet 197%================================================= 198% 199% 200% ** set up the conditional prior variance sg0bi and sgpbi. 201sg0bida = mu(1)^2*sg0bid; % ith equation 202sgpbida = mu(1)^2*mu(2)^2*sgpbid; 203sgpbida(ncoef-nexo+1:ncoef) = mu(1)^2*mu(3)^2; 204 %<<>> No scaling adjustment has been made for exogenous terms other than constant 205sgppbd = sgpbida(nvar+1:ncoef); % corresponding to A++, in a Sims-Zha paper 206 207Hptd = zeros(ncoef); 208Hptdi=Hptd; 209Hptd(ncoef,ncoef)=sgppbd(ncoef-nvar); 210Hptdinv(ncoef,ncoef)=1./sgppbd(ncoef-nvar); 211 % condtional on A0i, H_plus_tilde 212 213 214if nargin<13 % <<>>1 Default is no asymmetric information 215 asym0 = ones(nvar,nvar); % if not ones, then we have relative (asymmetric) tightness 216 asymp = ones(ncoef-1,nvar); % for A+. Column -- equation 217end 218 219%**** Asymmetric Information 220%asym0 = ones(nvar,nvar); % if not ones, then we have relative (asymmetric) tightness 221%asymp = ones(ncoef-1,nvar); % pp: plus without constant. Column -- equation 222%>>>>>> B: asymmetric prior variance for asymp <<<<<<<< 223% 224%for i = 1:lags 225% rowif = (i-1)*nvar+1; 226% rowil = i*nvar; 227% idmatw0 = 0.5; % weight assigned to idmat0 in the formation of asymp 228% if (i==1) 229% asymp(rowif:rowil,:)=(1-idmatw0)*ones(nvar)+idmatw0*idmat0; % first lag 230% % note: idmat1 is already transposed. Column -- equation 231% else 232% %asymp(rowif:rowil,1:nvar) = (1-idmatw0)*ones(nvar)+idmatw0*idmat0; 233% % <<<<<<< toggle + 234% % Note: already transposed, since idmat0 is transposed. 235% % Meaning: column implies equation 236% asymp(rowif:rowil,1:nvar) = ones(nvar); 237% % >>>>>>> toggle - 238% end 239%end 240% 241%>>>>>> E: asymmetric prior variance for asymp <<<<<<<< 242 243 244%================================================= 245% Computing the final covariance matrix (S1,...,Sm) for the prior of A0, 246% and final Gb=(G1,...,Gm) for A+ if asymmetric prior or for 247% B if symmetric prior for A+ 248%================================================= 249% 250for i = 1:nvar 251 %------------------------------ 252 % Introduce prior information on which variables "belong" in various equations. 253 % In this first trial, we just introduce this information here, in a model-specific way. 254 % Eventually this info has to be passed parametricly. In our first shot, we just damp down 255 % all coefficients except those on the diagonal. 256 257 %*** For A0 258 factor0=asym0(:,i); 259 260 sg0bd = sg0bida.*factor0; % Note, this only works for the prior variance Sg(i) 261 % of a0(i) being diagonal. If the prior variance Sg(i) is not 262 % diagonal, we have to the inverse to get inv(Sg(i)). 263 %sg0bdinv = 1./sg0bd; 264 % * unconditional variance on A0+ 265 H0td = diag(sg0bd); % unconditional 266 %=== Correlation in the MS equation to get a liquidity effect. 267 if (i==indxmsmdeqn(1)) 268 H0td(indxmsmdeqn(3),indxmsmdeqn(4)) = hpmsmd(1)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4))); 269 H0td(indxmsmdeqn(4),indxmsmdeqn(3)) = hpmsmd(1)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4))); 270 elseif (i==indxmsmdeqn(2)) 271 H0td(indxmsmdeqn(3),indxmsmdeqn(4)) = hpmsmd(2)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4))); 272 H0td(indxmsmdeqn(4),indxmsmdeqn(3)) = hpmsmd(2)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4))); 273 end 274 H0tdinv = inv(H0td); 275 %H0tdinv = diag(sg0bdinv); 276 % 277 H0tldcell_inv{i}=(Ui{i}'*kron(eye(nStates),H0tdinv/nStates))*Ui{i}; 278 279 280 %*** For A+ 281 if ~(lags==0) % For A1 to remain random walk properties 282 factor1=asymp(1:nvar,i); 283 sg1bd = sgpbida(1:nvar).*factor1; 284 sg1bdinv = 1./sg1bd; 285 % 286 Hptd(1:nvar,1:nvar)=diag(sg1bd); 287 Hptdinv(1:nvar,1:nvar)=diag(sg1bdinv); 288 if lags>1 289 factorpp=asymp(nvar+1:ncoef-1,i); 290 sgpp_cbd = sgppbd(1:ncoef-nvar-1) .* factorpp; 291 sgpp_cbdinv = 1./sgpp_cbd; 292 Hptd(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbd); 293 Hptdinv(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbdinv); 294 % condtional on A0i, H_plus_tilde 295 end 296 end 297 %--------------- 298 % The dummy observation prior affects only the prior covariance of A+|A0, 299 % but not the covariance of A0. See pp.69a-69b for the proof. 300 %--------------- 301 if indxDummy % Dummy observations as part of the explicit prior. 302 Hptdinv2 = Hptdinv + xtxbar; % Rename Hptdinv to Hptdinv2 because we want to keep Hptdinv diagonal in the next loop of i. 303 else 304 Hptdinv2 = Hptdinv; 305 end 306 Hptldcell_inv{i}=(Vi{i}'*kron(eye(nStates),Hptdinv2/nStates))*Vi{i}; 307 %Hptdinv_3 = kron(eye(nStates),Hptdinv); % ????? 308end 309 310 311