1function [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] ... 2 = fn_rnrprior_covres(nvar,q_m,lags,xdgel,mu,hpmsmd,indxmsmdeqn,nexo,asym0,asymp) 3% [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] ... 4% = fn_rnrprior_covres(nvar,q_m,lags,xdgel,mu,hpmsmd,indxmsmdeqn,nexo,asym0,asymp) 5% 6% More general than fn_rnrprior.m because when hpmsmd=0, fn_rnrprior_covres() is the same as fn_rnrprior(). 7% Allows for prior covariances for the MS and MD equations to achieve liquidity effects. 8% Exports random Bayesian prior of Sims and Zha with asymmetric rior (but no linear restrictions yet) 9% See Waggoner and Zha's Gibbs sampling paper and TVBVAR NOTES p. 71k.0. 10% 11% nvar: number of endogenous variables 12% q_m: quarter or month 13% lags: the maximum length of lag 14% xdgel: the general matrix of the original data (no manipulation involved) 15% with sample size including lags. Used only to get variances of residuals for 16% the scaling purpose; NOT used for mu(5) and mu(6). 17% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where 18% mu(5) and mu(6) are NOT used here. See fn_dataxy.m for using mu(5) and mu(6). 19% mu(1): overall tightness and also for A0; (0.57) 20% mu(2): relative tightness for A+; (0.13) 21% mu(3): relative tightness for the constant term; (0.1). NOTE: for other 22% exogenous terms, the variance of each exogenous term must be taken into 23% acount to eliminate the scaling factor. 24% mu(4): tightness on lag decay; (1) 25% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5) 26% mu(6): weight on single dummy initial observation including constant 27% (cointegration, unit roots, and stationarity); (5) 28% NOTE: for this subdirectory, mu(5) and mu(6) are not used. See fn_dataxy.m for using mu(5) and mu(6). 29% 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. 30% 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. 31% 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. 32% This will give us a liquidity effect. 33% indxmsmdeqn: 4-by-1 index for the locations of the MS and MD equation and for the locations of M and R. 34% indxmsmdeqn(1) for MS and indxmsmdeqn(2) for MD. 35% indxmsmdeqn(3) for M and indxmsmdeqn(4) for R. 36% nexo: number of exogenous variables (if not specified, nexo=1 (constant) by default). 37% The constant term is always put to the last of all endogenous and exogenous variables. 38% asym0: nvar-by-nvar asymmetric prior on A0. Column -- equation. 39% If ones(nvar,nvar), symmetric prior; if not, relative (asymmetric) tightness on A0. 40% asymp: ncoef-1-by-nvar asymmetric prior on A+ bar constant. Column -- equation. 41% If ones(ncoef-1,nvar), symmetric prior; if not, relative (asymmetric) tightness on A+. 42% -------------------- 43% Pi: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations 44% H0multi: nvar-by-nvar-by-nvar; H0 for different equations under asymmetric prior 45% Hpmulti: ncoef-by-ncoef-by-nvar; H+ for different equations under asymmetric prior 46% H0invmulti: nvar-by-nvar-by-nvar; inv(H0) for different equations under asymmetric prior 47% Hpinvmulti: ncoef-by-ncoef-by-nvar; inv(H+) for different equations under asymmetric prior 48% 49% Tao Zha, February 2000. Revised, September 2000, February 2003. 50% 51% Copyright (C) 1997-2012 Tao Zha 52% 53% This free software: you can redistribute it and/or modify 54% it under the terms of the GNU General Public License as published by 55% the Free Software Foundation, either version 3 of the License, or 56% (at your option) any later version. 57% 58% It is distributed in the hope that it will be useful, 59% but WITHOUT ANY WARRANTY; without even the implied warranty of 60% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 61% GNU General Public License for more details. 62% 63% If you did not received a copy of the GNU General Public License 64% with this software, see <http://www.gnu.org/licenses/>. 65% 66 67 68 69if nargin==7, nexo=1; end 70ncoef = nvar*lags+nexo; % number of coefficients in *each* equation, RHS coefficients only. 71 72H0multi=zeros(nvar,nvar,nvar); % H0 for different equations under asymmetric prior 73Hpmulti=zeros(ncoef,ncoef,nvar); % H+ for different equations under asymmetric prior 74H0invmulti=zeros(nvar,nvar,nvar); % inv(H0) for different equations under asymmetric prior 75Hpinvmulti=zeros(ncoef,ncoef,nvar); % inv(H+) for different equations under asymmetric prior 76 77%*** Constructing Pi for the ith equation under the random walk assumption 78Pi = zeros(ncoef,nvar); % same for all equations 79Pi(1:nvar,1:nvar) = eye(nvar); % random walk 80 81% 82%@@@ Prepared for Bayesian prior 83% 84% 85% ** monthly lag decay in order to match quarterly decay: a*exp(bl) where 86% ** l is the monthly lag. Suppose quarterly decay is 1/x where x=1,2,3,4. 87% ** Let the decay of l1 (a*exp(b*l1)) match that of x1 (say, beginning: 1/1) 88% ** and the decay of l2 (a*exp(b*l2)) match that of x2 (say, end: 1/5), 89% ** we can solve for a and b which are 90% ** b = (log_x1-log_x2)/(l1-l2), and a = x1*exp(-b*l1). 91if q_m==12 92 l1 = 1; % 1st month == 1st quarter 93 xx1 = 1; % 1st quarter 94 l2 = lags; % last month 95 xx2 = 1/((ceil(lags/3))^mu(4)); % last quarter 96 %xx2 = 1/6; % last quarter 97 % 3rd quarter: i.e., we intend to let decay of the 6th month match 98 % that of the 3rd quarter, so that the 6th month decays a little 99 % faster than the second quarter which is 1/2. 100 if lags==1 101 b = 0; 102 else 103 b = (log(xx1)-log(xx2))/(l1-l2); 104 end 105 a = xx1*exp(-b*l1); 106end 107 108 109 110% 111% *** specify the prior for each equation separately, SZ method, 112% ** get the residuals from univariate regressions. 113% 114sgh = zeros(nvar,1); % square root 115sgsh = sgh; % square 116nSample=size(xdgel,1); % sample size-lags 117yu = xdgel; 118C = ones(nSample,1); 119for k=1:nvar 120 [Bk,ek,junk1,junk2,junk3,junk4] = sye([yu(:,k) C],lags); 121 clear Bk junk1 junk2 junk3 junk4; 122 sgsh(k) = ek'*ek/(nSample-lags); 123 sgh(k) = sqrt(sgsh(k)); 124end 125% ** prior variance for A0(:,1), same for all equations!!! 126sg0bid = zeros(nvar,1); % Sigma0_bar diagonal only for the ith equation 127for j=1:nvar 128 sg0bid(j) = 1/sgsh(j); % sgsh = sigmai^2 129end 130% ** prior variance for lagged and exogeous variables, same for all equations 131sgpbid = zeros(ncoef,1); % Sigma_plus_bar, diagonal, for the ith equation 132for i = 1:lags 133 if (q_m==12) 134 lagdecay = a*exp(b*i*mu(4)); 135 end 136 % 137 for j = 1:nvar 138 if (q_m==12) 139 % exponential decay to match quarterly decay 140 sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation 141 elseif (q_m==4) 142 sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation 143 else 144 error('Incompatibility with lags, check the possible errors!!!') 145 %warning('Incompatibility with lags, check the possible errors!!!') 146 %return 147 end 148 end 149end 150% 151 152 153%================================================= 154% Computing the (prior) covariance matrix for the posterior of A0, no data yet 155%================================================= 156% 157% 158% ** set up the conditional prior variance sg0bi and sgpbi. 159sg0bida = mu(1)^2*sg0bid; % ith equation 160sgpbida = mu(1)^2*mu(2)^2*sgpbid; 161sgpbida(ncoef-nexo+1:ncoef) = mu(1)^2*mu(3)^2; 162 %<<>> No scaling adjustment has been made for exogenous terms other than constant 163sgppbd = sgpbida(nvar+1:ncoef); % corresponding to A++, in a Sims-Zha paper 164 165Hptd = zeros(ncoef); 166Hptdi=Hptd; 167Hptd(ncoef,ncoef)=sgppbd(ncoef-nvar); 168Hptdinv(ncoef,ncoef)=1./sgppbd(ncoef-nvar); 169 % condtional on A0i, H_plus_tilde 170 171 172if nargin<9 % the default is no asymmetric information 173 asym0 = ones(nvar,nvar); % if not ones, then we have relative (asymmetric) tightness 174 asymp = ones(ncoef-1,nvar); % for A+. Column -- equation 175end 176 177%**** Asymmetric Information 178%asym0 = ones(nvar,nvar); % if not ones, then we have relative (asymmetric) tightness 179%asymp = ones(ncoef-1,nvar); % pp: plus without constant. Column -- equation 180%>>>>>> B: asymmetric prior variance for asymp <<<<<<<< 181% 182%for i = 1:lags 183% rowif = (i-1)*nvar+1; 184% rowil = i*nvar; 185% idmatw0 = 0.5; % weight assigned to idmat0 in the formation of asymp 186% if (i==1) 187% asymp(rowif:rowil,:)=(1-idmatw0)*ones(nvar)+idmatw0*idmat0; % first lag 188% % note: idmat1 is already transposed. Column -- equation 189% else 190% %asymp(rowif:rowil,1:nvar) = (1-idmatw0)*ones(nvar)+idmatw0*idmat0; 191% % <<<<<<< toggle + 192% % Note: already transposed, since idmat0 is transposed. 193% % Meaning: column implies equation 194% asymp(rowif:rowil,1:nvar) = ones(nvar); 195% % >>>>>>> toggle - 196% end 197%end 198% 199%>>>>>> E: asymmetric prior variance for asymp <<<<<<<< 200 201 202%================================================= 203% Computing the final covariance matrix (S1,...,Sm) for the prior of A0, 204% and final Gb=(G1,...,Gm) for A+ if asymmetric prior or for 205% B if symmetric prior for A+ 206%================================================= 207% 208for i = 1:nvar 209 %------------------------------ 210 % Introduce prior information on which variables "belong" in various equations. 211 % In this first trial, we just introduce this information here, in a model-specific way. 212 % Eventually this info has to be passed parametricly. In our first shot, we just damp down 213 % all coefficients except those on the diagonal. 214 215 %*** For A0 216 factor0=asym0(:,i); 217 sg0bd = sg0bida.*factor0; % Note, this only works for the prior variance Sg(i) 218 % of a0(i) being diagonal. If the prior variance Sg(i) is not 219 % diagonal, we have to the inverse to get inv(Sg(i)). 220 %sg0bdinv = 1./sg0bd; 221 % * unconditional variance on A0+ 222 H0td = diag(sg0bd); % unconditional 223 %=== Correlation in the MS equation to get a liquidity effect. 224 if (i==indxmsmdeqn(1)) 225 H0td(indxmsmdeqn(3),indxmsmdeqn(4)) = hpmsmd(1)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4))); 226 H0td(indxmsmdeqn(4),indxmsmdeqn(3)) = hpmsmd(1)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4))); 227 else 228 H0td(indxmsmdeqn(3),indxmsmdeqn(4)) = hpmsmd(2)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4))); 229 H0td(indxmsmdeqn(4),indxmsmdeqn(3)) = hpmsmd(2)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4))); 230 end 231 H0tdinv = inv(H0td); 232 %H0tdinv = diag(sg0bdinv); 233 % 234 H0multi(:,:,i)=H0td; 235 H0invmulti(:,:,i)=H0tdinv; 236 237 238 %*** For A+ 239 if ~(lags==0) % For A1 to remain random walk properties 240 factor1=asymp(1:nvar,i); 241 sg1bd = sgpbida(1:nvar).*factor1; 242 sg1bdinv = 1./sg1bd; 243 % 244 Hptd(1:nvar,1:nvar)=diag(sg1bd); 245 Hptdinv(1:nvar,1:nvar)=diag(sg1bdinv); 246 if lags>1 247 factorpp=asymp(nvar+1:ncoef-1,i); 248 sgpp_cbd = sgppbd(1:ncoef-nvar-1) .* factorpp; 249 sgpp_cbdinv = 1./sgpp_cbd; 250 Hptd(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbd); 251 Hptdinv(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbdinv); 252 % condtional on A0i, H_plus_tilde 253 end 254 end 255 Hpmulti(:,:,i)=Hptd; 256 Hpinvmulti(:,:,i)=Hptdinv; 257end 258 259 260