1function of = pmddf236(x,idmat0,idmatp,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% Revsions by TZ, 10/13/96: efficiency improvement by streamlining the previous code 23% according to the general setup in Sims and Zha "Bayesian Methods for ...". 24% 25% Copyright (C) 1997-2012 Christopher A. Sims and Tao Zha 26% 27% This free software: you can redistribute it and/or modify 28% it under the terms of the GNU General Public License as published by 29% the Free Software Foundation, either version 3 of the License, or 30% (at your option) any later version. 31% 32% It is distributed in the hope that it will be useful, 33% but WITHOUT ANY WARRANTY; without even the implied warranty of 34% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 35% GNU General Public License for more details. 36% 37% If you did not received a copy of the GNU General Public License 38% with this software, see <http://www.gnu.org/licenses/>. 39% 40 41global vaha0 vaMMa0 GlbAph FRESHFUNCTION 42% FRESHFUNCTION added 8/20/90 by CAS. It signals the pmddg23 whether the global variables 43% it uses have been refreshed by a new pmddf231 call since the last gradient call. 44a0indx=find(idmat0); 45% pre-allocate GlbAph to be zeros(nvar,lags*nvar+1=ncoef) 46 47nhp = 0; % <<>> 4 hyperparameters 48na0p = length(a0indx); % <<>> number of A0 parameters 49nfp = na0p+nhp; 50 51% *** initializing 52vaha0 = zeros(na0p,1); 53vaMMa0 = zeros(na0p,1); 54 55% 56A0h = zeros(nvar,nvar); 57A0h(a0indx) = x(1:na0p); 58%A0h = reshape(a0,nvar,nvar); CAS 9/24/96. The reshape is not necessary if a0 starts as nvar x nvar 59 % restrictions on A0 60 61[A0hl,A0hu] = lu(A0h); 62 63% ** hyperparameters 64%%mu = x(na0p+1:nfp); 65mu = ones(5,1); 66mu(1) = 1; 67%mu(1) = 2; 68mu(2) = 0.2; 69mu(3) = 1; 70%mu(4)=1; 71mu(4)=10; 72mu(5) =1; 73%mu(5) = 40; 74% results from ...\dummy1\foreh\pmdd6.out 75% 76% mu(1): overall tightness and also for A0; 77% mu(2): relative tightness for A+; 78% mu(3): relative tightness for the constant term; 79% mu(4): weight on single dummy initial observation including constant; 80% mu(5): weight on nvar sums of coeffs dummy observations. 81 82% ** weight prior dummy observations 83%phi(1:nvar,:) = (mu(5)^2/mu(1)^2)*phi(1:nvar,:); 84%y(1:nvar,:) = (mu(5)^2/mu(1)^2)*y(1:nvar,:); 85%phi(ndobs,:) = mu 86% modified by CAS 8/6/96. The dummy obs weights are naturally in std units, not var units. 87phi(1:nvar,:) = mu(5)*phi(1:nvar,:); 88y(1:nvar,:) = mu(5)*y(1:nvar,:); 89phi(nvar+1,:) = mu(4)*phi(nvar+1,:); 90y(nvar+1,:) = mu(4)*y(nvar+1,:); 91 92% ** set up the conditional prior variance sg0bi and sgpbi. 93sg0bida = mu(1)^2*sg0bid; 94sgpbida = mu(1)^2*mu(2)^2*sgpbid; 95sgpbida(ncoef) = mu(1)^2*mu(3)^2; 96sgppbd = sgpbida(nvar+1:ncoef); % corresponding to A++, in an SZ paper 97%%sgppbdi = 1./sgppbd; 98%sgppb = 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 109yty=y'*y; 110ymy=yty-cxyxx; %trial 111%ymy = y'*M*y; 112%cxy = phi'*y; %trial 113cxy = vd*yu'; %trial 114cyx = cxy'; 115%cxpy = xtx\cxy; %not used further except in line below 116%cxyxx = cxy'*cxpy; 117 118% ** estimate A+_hat conditional on A0, ncoef*nvar, but not using full prior 119%GlbAph = cxpy*A0h; 120 121 122%%%%%%% 123%%% build the objective function below 124%%%%%%% 125% 126t1f = diag(abs(A0hu)); 127t1f = log(t1f); 128% * without the prior |A0|^k 129tt1 = (-1.0)*fss*sum(t1f); 130 131tt3 = 0; 132 133disp(sprintf('Starting loop (of): %g',toc)) 134 135stril = 0; 136Hptd = zeros(ncoef); 137Hptdi=Hptd; 138%%Hptd(nvar+1:ncoef,nvar+1:ncoef)=diag(sgppbd); 139%%Hptdi(nvar+1:ncoef,nvar+1:ncoef)=diag(sgppbdi); 140Hptd(ncoef,ncoef)=sgppbd(ncoef-nvar); 141Hptdi(ncoef,ncoef)=1/sgppbd(ncoef-nvar); 142 % condtional on A0i, H_plus_tilde 143for i = 1:nvar 144 stri = find(idmat0(:,i)); %CAS change 9/24/96, to avoid slim chance of an exact zero in x vector 145 %at some point during an optimization. 146 strm = length(stri); 147 148 A0hb = A0h(stri,i)-A0b(stri,i); 149 150 % ** set up the conditional prior variance sg0bi and sgpbi. 151 % sg0bd = zeros(stri,1); Commented out by CAS, 9/24/96. Looks like this meant to have 152 % strm where it has stri, and in any case it is "overwritten" below. 153 %------------------------------ 154 % Introduce prior information on which variables "belong" in various equations. 155 % In this first trial, we just introduce this information here, in a model-specific way. 156 % Eventually this info has to be passed parametricly. In our first shot, we just damp down 157 % all coefficients except those on the diagonal. 158 factor0=idmat0(:,i); 159 sg0bd = sg0bida(stri).*factor0(stri); 160 sg0bdi = 1./sg0bd; 161 %sg0b = diag(sg0bd); 162 % 163 factor1=idmatp(1:nvar,i); 164 sg1bd = sgpbida(1:nvar).*factor1; 165 sg1bdi = 1./sg1bd; 166 %sg1b = diag(sg1bd); 167 % 168 factorpp=idmatp(nvar+1:ncoef-1,i); 169 sgpp_cbd = sgppbd(1:ncoef-nvar-1) .* factorpp; 170 sgpp_cbdi = 1./sgpp_cbd; 171 172 173 % ** set up the unconditional prior variance on A0i and A+i 174 XX = zeros(nvar,strm); 175 XX(stri,:) = eye(strm); % XX in Gelman, etel, on p479 176 % 177 % * final conditional prior variance on A+i give that on A0i, and 178 % * unconditional variance on A0+ 179 H0td = diag(sg0bd); % unconditional 180 % H_~: td: tilde for ~ 181 % ** inverse and chol decomp 182 H0tdi = diag(sg0bdi); 183 184 % 185 Hptd(1:nvar,1:nvar)=diag(sg1bd); 186 Hptdi(1:nvar,1:nvar)=diag(sg1bdi); 187 Hptd(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbd); 188 Hptdi(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbdi); 189 % condtional on A0i, H_plus_tilde 190 191 % common terms 192 xxhp = xtx+Hptdi; 193 A0hbH = A0hb'*H0tdi; 194 H1p_1 = zeros(ncoef,nvar); 195 H1p_1(1:nvar,:) = Hptdi(1:nvar,1:nvar); 196 %%Hm = (cyx+H1p_1')*(xxhp\(cxy+H1p_1)); 197 Hm1 = (cyx+H1p_1')/xxhp; 198 Hm2 = cxy+H1p_1; 199 GlbAph(i,:) = A0h(stri,i)'*XX'*Hm1; % alpha_plus_*_transpose 200 %%alpMpyh = A0h(stri,i)'*XX'*(yty+Hptdi(1:nvar,1:nvar)-Hm)*XX; 201 alpMpyh = A0h(stri,i)'*XX'*(yty+Hptdi(1:nvar,1:nvar))*XX - GlbAph(i,:)*Hm2*XX; 202 203 % ** 3rd bid term in i-th equation, 204 % ** with all other 3rd big terms before the i-th equation 205 tt3 = tt3 + A0hbH*A0hb + alpMpyh*A0h(stri,i); 206 207 %%%% 208 % *** passing the gradient to "pmdg6.m" (analytical gradient) 209 %%%% 210 % 211 % ** daha0 = d((alpha0-alpha0_tilde)'*inv(H_0_tilde)* 212 % ** (alpha0-alpha0_tilde))/d(a0??) 213 daha0 = 2.0*A0hbH; 214 vaha0(stril+1:stril+strm) = daha0(:); 215 % 216 % ** daMMa0 = d(alpha0'*P'*(Y'Y+H_1^(-1)+H_m)*P*apha0) / d(a0?) 217 daMMa0 = 2.0*alpMpyh; 218 vaMMa0(stril+1:stril+strm) = daMMa0(:); 219 220 stril = stril + strm; 221end 222 223%disp(sprintf('Loop %d end: %g', i, toc)) 224disp(sprintf('Loop end (of): %g', toc)) 225 226%e_tfat = toc 227 228of = tt1 + 0.5*tt3; 229FRESHFUNCTION=1; 230