1function of = pmddf233(x,idmat0,idmat1,fss,nvar,ncoef,phi, ... 2 y,A0b,sg0bid,sgpbid) 3%function of = pmddf232(x,idmat,fss,nvar,ncoef,phi, ... 4% y,A0b,sg0bid,sgpbid) 5% Leeper-Sims-Zha BVAR setup 6% general program to setup A0 matrix and compute the likelihood 7% requires 8% x (parameter vector) 9% a0indx (matrix indicating the free parameters in A0) 10% fss (forecast sample size) 11% nvar (number of variables) 12% ncoef (number of coefficients in a single equation in A+) 13% phi (r.h.s. observations in the system for A+) 14% y (l.h.s. observations for A0) 15% A0b (initial prior mean on A0) 16% sg0bid (diagonal of Sigma0_bar, unweighted prior vcv on the parameters in i-th equation in A0) 17% sgpbid (diagonal of Sigma+_bar, unweighted prior vcv on the parameters in i-th equation in A+) 18% 19% Revisions by CAS 9/27/96: If we transmit idmat, rather than a0indx, we can scale the elements of 20% idmat so that it carries information about relative prior variances. 21 22% 23% Copyright (C) 1997-2012 Christopher A. Sims 24% 25% This free software: you can redistribute it and/or modify 26% it under the terms of the GNU General Public License as published by 27% the Free Software Foundation, either version 3 of the License, or 28% (at your option) any later version. 29% 30% It is distributed in the hope that it will be useful, 31% but WITHOUT ANY WARRANTY; without even the implied warranty of 32% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 33% GNU General Public License for more details. 34% 35% If you did not received a copy of the GNU General Public License 36% with this software, see <http://www.gnu.org/licenses/>. 37% 38 39global vaya0 vaha0 vaxah vahad vbiga GlbAph FRESHFUNCTION 40% FRESHFUNCTION added 8/20/90 by CAS. It signals the pmddg23 whether the global variables 41% it uses have been refreshed by a new pmddf231 call since the last gradient call. 42a0indx=find(idmat0); 43% pre-allocate GlbAph to be zeros(nvar,lags*nvar+1=ncoef) 44 45nhp = 0; % <<>> 4 hyperparameters 46na0p = length(a0indx); % <<>> number of A0 parameters 47nfp = na0p+nhp; 48 49% *** initializing 50vaya0 = zeros(na0p,1); 51vaha0 = zeros(na0p,1); 52vaxah = zeros(na0p,1); 53vahad = zeros(na0p,1); 54vbiga = zeros(na0p,1); 55 56 57A0h = zeros(nvar,nvar); 58A0h(a0indx) = x(1:na0p); 59%A0h = reshape(a0,nvar,nvar); CAS 9/24/96. The reshape is not necessary if a0 starts as nvar x nvar 60 % restrictions on A0 61 62[A0hl,A0hu] = lu(A0h); 63 64% ** hyperparameters 65%%mu = x(na0p+1:nfp); 66mu = ones(5,1); 67mu(1) = 1; 68%mu(1) = 2; 69mu(2) = 0.2; 70mu(3) = 1; 71%mu(4)=1; 72mu(4)=1; 73mu(5) =1; 74%mu(5) = 40; 75% results from ...\dummy1\foreh\pmdd6.out 76% 77% mu(1): overall tightness and also for A0; 78% mu(2): relative tightness for A+; 79% mu(3): relative tightness for the constant term; 80% mu(4): weight on single dummy initial observation including constant; 81% mu(5): weight on nvar sums of coeffs dummy observations. 82 83% ** weight prior dummy observations 84%phi(1:nvar,:) = (mu(5)^2/mu(1)^2)*phi(1:nvar,:); 85%y(1:nvar,:) = (mu(5)^2/mu(1)^2)*y(1:nvar,:); 86%phi(ndobs,:) = mu 87% modified by CAS 8/6/96. The dummy obs weights are naturally in std units, not var units. 88phi(1:nvar,:) = mu(5)*phi(1:nvar,:); 89y(1:nvar,:) = mu(5)*y(1:nvar,:); 90phi(nvar+1,:) = mu(4)*phi(nvar+1,:); 91y(nvar+1,:) = mu(4)*y(nvar+1,:); 92 93% ** set up the conditional prior variance sg0bi and sgpbi. 94sg0bida = mu(1)^2*sg0bid; 95sgpbida = mu(1)^2*mu(2)^2*sgpbid; 96sgpbida(ncoef) = mu(1)^2*mu(3)^2; 97sgppbd = sgpbida(nvar+1:ncoef); % corresponding to A++, in a SZ paper 98sgppb = diag(sgppbd); 99 100% 101% ** some common terms 102[u d v]=svd(phi,0); %trial 103% xtx = phi'*phi; %trial 104vd=v.*(ones(size(v,2),1)*diag(d)'); %trial 105xtx=vd*vd'; 106% M = eye(fss) - phi*(xtx\phi'); % not used except in line below %trial 107yu = y'*u; %trial 108cxyxx=yu*yu'; %trial 109ymy=y'*y-cxyxx; %trial 110%ymy = y'*M*y; 111%cxy = phi'*y; %trial 112cxy = vd*yu'; %trial 113%cxpy = xtx\cxy; %not used further except in line below 114%cxyxx = cxy'*cxpy; 115 116% ** estimate A+_hat conditional on A0, ncoef*nvar, but not using full prior 117%GlbAph = cxpy*A0h; 118 119 120%%%%%%% 121%%% build the objective function below 122%%%%%%% 123% 124t1f = diag(abs(A0hu)); 125t1f = log(t1f); 126tt1 = (-1.0)*fss*sum(t1f); 127% Commented out by TZ, 10/3/96, to allow the prior |A0|^k 128%%tt1 = (-1.0)*(fss+ncoef)*sum(t1f); 129 130tt3 = 0; 131tt4 = 0; 132 133disp(sprintf('Starting loop (of): %g',toc)) 134 135stril = 0; 136Hp0p = zeros(ncoef); 137% initializing Htd_+0*inv(Htd_00)*Htd_0+: conditional prior covariance 138for i = 1:nvar 139 stri = find(idmat0(:,i)); %CAS change 9/24/96, to avoid slim chance of an exact zero in x vector 140 %at some point during an optimization. 141 strm = length(stri); 142 143 A0hb = A0h(stri,i)-A0b(stri,i); 144 Hp0 = zeros(ncoef,strm); 145 % initializing Htd_+0*inv(Htd_00)*Htd_0 146 147 % ** set up the conditional prior variance sg0bi and sgpbi. 148 % sg0bd = zeros(stri,1); Commented out by CAS, 9/24/96. Looks like this meant to have 149 % strm where it has stri, and in any case it is "overwritten" below. 150 %------------------------------ 151 % Introduce prior information on which variables "belong" in various equations. 152 % In this first trial, we just introduce this information here, in a model-specific way. 153 % Eventually this info has to be passed parametricly. In our first shot, we just damp down 154 % all coefficients except those on the diagonal. 155 factor0=idmat0(:,i); 156 sg0bd = sg0bida(stri).*factor0(stri); 157 sg0b = diag(sg0bd); 158 % 159 factor1=idmat1(:,i); 160 sg1bd = sgpbida(1:nvar).*factor1; 161 sg1b = diag(sg1bd); 162 163 % ** set up the unconditional prior variance on A0i and A+i 164 XX = zeros(nvar,strm); 165 XX(stri,:) = eye(strm); % XX in Gelman, etel, on p479 166 % 167 Hb = zeros(strm+ncoef); % including A0i and A+i 168 Hb(1:strm,1:strm) = sg0b; 169 sg0bxx = sg0b*XX'; 170 Hb(1:strm,strm+1:strm+nvar) = sg0bxx; 171 %Hb(strm+1:strm+nvar,1:strm) = XX*sg0b; CAS 9/24/96. Saves a little multiplication by using line below. 172 Hb(strm+1:strm+nvar,1:strm) = sg0bxx'; 173 Hb(strm+1:strm+nvar,strm+1:strm+nvar) = XX*sg0bxx+sg1b; 174 Hb(strm+nvar+1:strm+ncoef,strm+nvar+1:strm+ncoef) = sgppb; 175 176 % ** set up the final prior variance on A0i and A+i, i =1,..,6 separately. 177 % ** using dummy initial observations 178 %zs1 = [ys1(:,stri) -xd0]; % z*1 for the 1st equation 179 180 % * final unconditional prior variance on A0i and A+i 181 %Hbzs1=Hb*zs1'; 182 %Htd = Hb - Hbzs1*((zs1*Hbzs1+mu(1)^2*mu(4)^2*eye(nvar))\Hbzs1'); 183 % CAS mod 8/7/96: Hope is that dropping above two lines in favor of the one below 184 % removes double-counting of dummy observations and a possible scale sensitivity. 185 Htd = Hb; 186 % H_~: td: tilde for ~ 187 188 % * final conditional prior variance on A+i give that on A0i, and 189 % * unconditional variance on A0+ 190 H0td = Htd(1:strm,1:strm); % unconditional 191 % ** inverse and chol decomp 192 H0tdi = inv(H0td); 193 194 % 195 Hptd10 = Htd(strm+1:strm+nvar,1:strm); 196 % top matrix of nvar in Htd_+0 197 Hptd100 = Hptd10*H0tdi; 198 Hp0(1:nvar,:) = Hptd100; 199 Hptd101 = Hptd100*Hptd10'; 200 Hp0p(1:nvar,1:nvar) = Hptd101; 201 Hptd = Htd(strm+1:strm+ncoef,strm+1:strm+ncoef) - Hp0p; 202 % condtional on A0i, H_plus_tilde 203 204 % ** inverse 205 Hptdi = inv(Hptd); 206 chhh = Hptdi*Hp0; % <<>> the term not used for "of" 207 208 % common term 209 xxhp = xtx+Hptdi; 210 211 % ** final conditional prior mean on A+i 212 alptd = zeros(ncoef,1); 213 alptd(1:nvar) = XX*A0b(stri,i); 214 %alptd = alptd + Hptdf*(H0tdi*A0hb); CAS efficiency improvement 8/8/96 215 alptd = alptd + Hp0*A0hb; 216 % conditional on A0i, alpha-plus-tilde 217 218 % * alpha_plus_* for the ith equation, ncoef*1 219 % *alps = xxhp\(phi'*y*A0h(:,i)+HptdDinalptd); 220 HptdDinalptd=Hptdi*alptd; 221 xaha = cxy*A0h(:,i)+HptdDinalptd; %CAS 8/7/96. Replaced phi'*y with cxy 222 alpst = xaha'/xxhp; 223 GlbAph(i,:)=alpst; 224 % alpst': transpose of alps. 225 226 % ** 3rd bid term in i-th equation, 227 % ** with all other 3rd big terms prior to i-th 228 A0hy = A0h(:,i)'*ymy; 229 A0hbH = A0hb'*H0tdi; 230 tt3 = tt3 + A0hy*A0h(:,i) + A0hbH*A0hb; 231 % ** 4th bid term in i-th equation, 232 % ** with all other 4th big terms prior to i-th 233 A0hcxyxx = A0h(:,i)'*cxyxx; 234 atdHatd = alptd'*HptdDinalptd; 235 tt4 = tt4 + A0hcxyxx*A0h(:,i)+atdHatd-alpst*xaha; 236 237 %%%% 238 % *** passing the gradient to "pmdg6.m" (analytical gradient) 239 %%%% 240 % ** d1aya0 = d(alpha0'*Y'MY*alpha0)/d(a0??) 241 daya0 = 2.0*A0hy; 242 daya0 = daya0(:); 243 vaya0(stril+1:stril+strm) = daya0(stri); 244 245 % ** daha0 = d((alpha0-alpha0_tilde)'*inv(H_0_tilde)* 246 % ** (alpha0-alpha0_tilde))/d(a0??) 247 daha0 = 2.0*A0hbH; 248 vaha0(stril+1:stril+strm) = daha0(:); 249 250 % ** daxah = d(alpha_+^'*X'X*alpha_+^)/d(a0??) 251 % ** daxah = d(alpha_+^'*X'X*alpha_+^)/d(a0??) 252 %% daxah = 2*(Aph(:,i)'*xtx) * cxpy(:,stri); 253 %% vaxah(stril+1:stril+strm) = daxah'; 254 daxah = 2*A0hcxyxx; 255 daxah = daxah(:); 256 vaxah(stril+1:stril+strm) = daxah(stri); 257 258 259 % ** dahad = d(alpha_+~'*inv(H_+~)*alpha_+~)/d(a0??) 260 dahad = 2*HptdDinalptd' * Hp0; 261 % <<>> multiplications not used in "of" 262 vahad(stril+1:stril+strm) = dahad'; 263 264 % ** dbiga = d((X'X*alpha_+^+inv(H_+~)*alpha_+*)'*a_+_*)/d(a0??) 265 dbiga = 2*alpst*(cxy(:,stri)+chhh); 266 % <<>> multiplications not used in "of" 267 vbiga(stril+1:stril+strm) = dbiga'; 268 269 stril = stril + strm; 270end 271 272%disp(sprintf('Loop %d end: %g', i, toc)) 273disp(sprintf('Loop end (of): %g', toc)) 274 275%e_tfat = toc 276 277of = tt1 + 0.5*tt3 + 0.5*tt4; 278FRESHFUNCTION=1; 279