1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
|
function imfsim(xinput)
% imfsim(xinput)
% Save the simulated pdfs of impulse responses;
% Print and save Gelman's measures of B and W for A0's,
% only when nstarts (# of starting points) >1.
% Ref: Waggoner and Zha "Does Normalization Matter for Inference"
% See note Forecast (2)
%
% xinput{1}: nfp -- total number of free parameters
% xinput{2}: nvar -- number of variables
% xinput{3}: xhat -- ML estimate of free parameters in A0
% xinput{4}: hess -- Hessian of -logLH
% xinput{5}:Indxv -- index for selected variables of interest; normall first 2 are of our interest
% to select variables, always check idmat0 to make sure
% it plots: (1) pdf of 1st v for every buffer, (2) scattered plot of 1st and 2nd for every buffer,
% (2) pdf of 1st v for all sequences; (4) scattered plot of 3rd and 4th for all sequences
% (5) scattered plot of 1st and 2nd for al sequences.
% xinput{6}: IndxGraph - 1: plot graphs; 0: no graphs
% xinput{7}: idmat0 -- Index for non-zero elements in A0 with column to equation
% xinput{8}: nstarts -- # of starting points in Gibbs sampling
% xinput{9}: ndraws1 -- # of 1st loop to be discarded
% xinput{10}: ndraws2 -- # of 2nd loop for saving A0 draws
% xinput{11}: imndraws=nstarts*ndraws2
% xinput{12}: a0indx -- index number for non-zero elements in A0
% xinput{13}: tdf -- degrees of freedom for t-distribution
% xinput{14}: nbuffer -- interval for printing, plotting, and saving
% xinput{15}: Sbd -- nvar-by-nvar S{1}, ..., S{n} -- kind of covariance matrix for each simultaneous equation
% Already divided by "fss."
% xinput{16}: nSample -- the original sample size including lags
% xinput{17}: IndxNmlr -- index for which normalization rule to choose
% xinput{18}: IndxGibbs -- index for WZ Gibbs; 1; Gibbs; 0: Metropolis
% xinput{19}: scf -- reduction scale factor for Metropolis jumping kernel
% xinput{20}: H_sr -- square root of the inverse of the covariance matrix
% for free elements in A0 (nfp-by-nfp)
% xinput{21}: fss -- effective sample size == nSample-lags+# of dummy observations
% xinput{22}: idfile1 -- calls "iden6std." Save stds. of both data and impulse responses in idfile1
% xinput{23}: xxhpc -- chol(X'X+inv(H_p_tilde)): upper triangular but its transpose
% is lower triagular Choleski
% xinput{24}: ImfErr -- if 1, impulse response simulation; if 0, disable this simulation
% xinput{25}: ninv -- number of bins pre-specified to put each draw of impulse response
% into a proper bin (or small interval)
% xinput{26}: imstp -- # of steps for impulse responses
% xinput{27}: forep -- forecast periods (# of steps)
% xinput{28}: yact -- actual data (in log except R, U, etc.)
% xinput{29}: yactqg -- quarterly annualized growth in actual data
% xinput{30}: yactCalyg -- calendar annual growth in actual data
% xinput{31}: imfml -- imstp-by-nvar^2 ML impulse responses
% xinput{32}: forepq -- forecast periods for quarterly growth
% xinput{33}: forepy -- forecast periods for annual growth
% xinput{34}: ncoef -- k: # of coeffients per equation
% xinput{35}: Bhml -- ML reduced form parameter B (nvar-by-k)
% xinput{36}: lags -- # of lags
%------------------
% imfcntmulti: (ninv+2,imstp*nvar^2,nstarts) matrix
% cnt: count; for impulse responses; multi (nstarts) sequences
% All output is saved in outB_W, including Range5, invc, ninv, and imfcntmulti
%
% Written by T. Zha 1999
%
% Copyright (C) 1997-2012 Tao Zha
%
% This free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% It is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% If you did not received a copy of the GNU General Public License
% with this software, see <http://www.gnu.org/licenses/>.
%
nfp=xinput{1}; nvar=xinput{2}; xhat=xinput{3}; hess=xinput{4}; Indxv=xinput{5};
IndxGraph=xinput{6}; idmat0=xinput{7}; nstarts=xinput{8}; ndraws1=xinput{9}; ndraws2=xinput{10};
imndraws=xinput{11}; a0indx=xinput{12}; tdf=xinput{13}; nbuffer=xinput{14}; Sbd=xinput{15};
nSample = xinput{16}; IndxNmlr=xinput{17}; IndxGibbs=xinput{18}; scf=xinput{19}; H_sr=xinput{20};
fss=xinput{21}; idfile1=xinput{22}; xxhpc=xinput{23}; ImfErr=xinput{24}; ninv=xinput{25};
imstp=xinput{26}; forep=xinput{27}; yact=xinput{28}; yactqg=xinput{29}; yactCalyg=xinput{30};
imfml=xinput{31}; forepq=xinput{32}; forepy=xinput{33}; ncoef=xinput{34}; Bhml=xinput{35};
lags=xinput{36};
Avhx_bs = zeros(nfp,1);
Avhx_bm = zeros(nfp,1);
Avhx_bj = zeros(nfp,1);
Avhx_cj = zeros(nfp,1);
Avhx_cm = zeros(nfp,1);
A0_h = zeros(nvar);
A0gbs = A0_h; % drawn A0 in Gibbs sampling
Avhxm = zeros(nfp,1);
Avhxs = Avhxm;
A0xhat = zeros(nvar);
A0xhat(a0indx) = xhat;
% A0hatw = zeros(nvar^2,nbuffer);
countJump = zeros(nstarts,1);
imfmean = zeros(imstp,nvar^2);
imfcntmulti = zeros(ninv+2,imstp*nvar^2,nstarts);
% cnt: count; for impulse responses; multi (nstarts) sequences
%---------------------------------------------------
% Specify the range for counting the empirical distribution
%
%** load the standard deviations of 6 variables, one for log(y), one for gq, and
%** the third one for yg
eval(['load ' idfile1 '.prn -ascii']);
eval(['ABstd=' idfile1 ';']);
Range5 = cell(4,1); % 4: log, qg, yg, and imf
%@@@ Tony's trick to expand the matrix
%
%** In order of log(y), qg, and yg for Range5{i} for i=1:3
Range5{1} =zeros(forep,nvar,2); % 2: min and max
Range5{1}(:,:,1) = repmat(yact(length(yact(:,1)),:)-10*ABstd(1,:),[forep 1]); % min, 30 std.
Range5{1}(:,:,2) = repmat(yact(length(yact(:,1)),:)+10*ABstd(1,:),[forep 1]); % max, 30 std.
%
Range5{2} =zeros(forepq,nvar,2); % 2: min and max
Range5{2}(:,:,1) = repmat(yactqg(length(yactqg(:,1)),:)-10*ABstd(2,:),[forepq 1]); % min, 30 std.
Range5{2}(:,:,2) = repmat(yactqg(length(yactqg(:,1)),:)+10*ABstd(2,:),[forepq 1]); % max, 30 std.
%
Range5{3} =zeros(forepy,nvar,2); % 2: min and max
Range5{3}(:,:,1) = repmat(yactCalyg(length(yactCalyg(:,1)),:)-10*ABstd(3,:),[forepy 1]); % min, 30 std.
Range5{3}(:,:,2) = repmat(yactCalyg(length(yactCalyg(:,1)),:)+10*ABstd(3,:),[forepy 1]); % max, 30 std.
%
Range5{4} =zeros(imstp,nvar^2,2); % 2: min and max
imfscale = repmat(ABstd(4,:),[1 nvar]); % because nvar variables to 1, ..., nvar shocks
Range5{4}(:,:,1) = repmat(imfml(1,:)-5*imfscale,[imstp 1]); % min, 5 std.
Range5{4}(:,:,2) = repmat(imfml(1,:)+5*imfscale,[imstp 1]); % max, 5 std.
% Range5(4)(:,:,1): imstp-by-nvar^2. Column: nvar responses to 1st shock,
% nvar responses to 2nd shock, ...
%**
invc = cell(4,1); % interval length (used for counting later). 1st 3 cells have each
% forep-by-nvar, and 4th has imstp-by-nvar^2.
for i=1:4
invc{i} = Range5{i}(:,:,2) - Range5{i}(:,:,1);
end
hbin = invc{4} ./ ninv; % bin size for each point of impulse responses kkdf
imfloor = Range5{4}(:,:,1);
%===================================
% Here begins with the big loop
%===================================
H1 = chol(hess); % upper triangular so that H1' is a lower triangular decomp
baseW = H_sr; %inv(H1); %H_sr; % covariance matrix without scaling
nswitch=0; %<<>> total number of sign switches
A0inxhat = inv(A0xhat); % inverse of ML estimate
a0indx0 = find(idmat0==0); % index for all zero's in A0;
nn=[nvar lags imstp];
[cT,vR,kdf] = gibbsglb(Sbd,idmat0,nvar,fss);
tic
for starts = 1:nstarts
starts
if starts == 1
A0gbs(a0indx) = xhat; % from "load ..."
if ~IndxGibbs % Metropolist
Avhx = xhat;
hAvhx = a0asfun(Avhx,Sbd,fss,nvar,a0indx);
hAvhx = -hAvhx; % converted to logLH
end
else
Avhx = baseW*randn(nfp,1); %H_sr*randn(nfp,1); % D: discarded sequence
csq=randn(tdf,1);
csq=sum(csq .* csq);
Avhx = xhat+Avhx/sqrt(csq/tdf);
%** Normalization by the choice of IndxNmlr
A0gbs(a0indx) = Avhx;
if ~IndxNmlr(5)
[A0gbs,jnk] = nmlzvar(A0gbs,A0xhat,A0inxhat,IndxNmlr,nswitch,[]);
else
A0ingbs = inv(A0gbs);
[A0gbs,jnk,jnk1] = nmlzvar(A0gbs,A0xhat,A0inxhat,IndxNmlr,nswitch,A0ingbs);
end
%
if ~IndxGibbs % Metropolist
Avhx = A0gbs(a0indx);
hAvhx = a0asfun(Avhx,Sbd,fss,nvar,a0indx);
hAvhx = -hAvhx; % converted to logLH
end
end
%
Avhxmm = zeros(nfp,1);
Avhxss = zeros(nfp,1);
cJump = 0;
imfcnt = zeros(ninv+2,imstp*nvar^2); % cnt: count; for impulse responses
for draws = 1:ndraws1
if IndxGibbs
A0gbs = gibbsvar(A0gbs,cT,vR,nvar,fss,kdf);
else % Metropolis
[Avhx,hAvhx,cJump] = smtplis(Avhx,hAvhx,tdf,cJump,scf,...
baseW,nfp,Sbd,fss,nvar,a0indx);
end
end
wdraws=(starts-1)*ndraws2+0;
for draws = 1:ndraws2
drawsc = (starts-1)*ndraws2+draws;
if IndxGibbs
A0gbs = gibbsvar(A0gbs,cT,vR,nvar,fss,kdf);
A0gbs(a0indx0) = 0; % set all zeros in A0gbs clean to avoid possible cumulative round-off errors
else % Metropolis
[Avhx,hAvhx,cJump] = smtplis(Avhx,hAvhx,tdf,cJump,scf,...
baseW,nfp,Sbd,fss,nvar,a0indx);
A0gbs(a0indx) = Avhx;
end
%*** call normalization so that A0_h is normalized
if ~IndxNmlr(5)
[A0_h,nswitch] = nmlzvar(A0gbs,A0xhat,A0inxhat,IndxNmlr,nswitch,[]);
A0_hin = inv(A0_h);
else
A0ingbs = inv(A0gbs);
[A0_h,nswitch,A0_hin] = nmlzvar(A0gbs,A0xhat,A0inxhat,IndxNmlr,nswitch,A0ingbs);
end
Avhx_norm = A0_h(a0indx);
%
% *** normal draws for posterior Aplus conditional on A0h
%
%** the mean is Aplushm, and the covariance is inv(xxhp)=Lxxhpc*Lxxhpc'
Apindm = randn(ncoef,nvar);
%
if ~all(all(finite(Bhml)))
Aplushm=zeros(ncoef,nvar);
for i=1:nvar
Aplushm(:,i)=Gb{i}*A0_h(:,i); % see Zha's forecast (1) p.9
% Here, Gb is used to compute A+ where A+(i) = Gb(i)*a0(i)
end
Bh_h = (Aplushm + xxhpc\Apindm)*A0_hin;
else
Bh_h = Bhml + (xxhpc\Apindm)*A0_hin;
end
if ImfErr
swish_h = A0_hin'; % Switching back to the form A0*y(t) = e(t)
imf_h = zimpulse(Bh_h,swish_h,nn); % in the form that is congenial to RATS
% imf3_h=reshape(imf_h,size(imf_h,1),nvar,nvar);
% % imf3: row--steps, column--nvar responses, 3rd dimension--nvar shocks
imfmean = imfmean + imf_h; % posterior mean
imfcnt = empdfsort(imfcnt,imf_h,imfloor,hbin,ninv);
% sorted counts (prob.) in bins
end
Avhxm = Avhxm + Avhx_norm; % 1st step to overall mean of parameter
Avhxs = Avhxs + Avhx_norm.^2; % 1st step to overall 2nd moment of parameter
%* compute the mean and 2nd moment
%** Getting average of variances W and variance of means B/n -- B_n
%* see Gelman, p.331, my Shock(0), 12-13, and my Forecast (2), 28-31
if (nstarts>1)
Avhxmm = Avhxmm + Avhx_norm; % n*(phi_.j)
Avhxss = Avhxss + Avhx_norm.^2; % 1st step to (phi_ij) for fixed j
end
% A0hatw(:,drawsc-wdraws) = A0_h(:);
if ~mod(draws,nbuffer)
starts
draws
wdraws=drawsc
% fwriteid = fopen('outA0.bin','a');
% count = fwrite(fwriteid,A0hatw,'double');
% status = fclose('all');
end
end
%
imfcntmulti(:,:,starts) = imfcnt;
if ~IndxGibbs
countJump(starts,1) = cJump;
end
%
%** Getting average of variances W and variance of means B/n -- B_n
%** see Gelman, p.331, my Shock(0), pp.12-13, and my Forecast (2), pp.28-31
if (nstarts>1)
Avhx_aj = Avhxmm/ndraws2; % (phi_.j)
Avhx_bs = Avhx_bs + Avhx_aj.^2; % 1st step to 2nd moment of (phi_.j)
Avhx_bm = Avhx_bm + Avhx_aj; % 1st step to (phi_..)
%
Avhx_bj = Avhxss/ndraws2; % 2nd moment of (phi_ij) for fixed j
Avhx_cj = Avhx_bj - Avhx_aj.^2; % ((n-1)/n)*S^2_j
Avhx_cm = Avhx_cm + Avhx_cj; % sum of ((n-1)/n)*S^2_j across j
end
end
timend = toc
timeminutes=timend/60
if ~IndxGibbs
countJump = countJump/ndraws2
end
Avhxm = Avhxm/(imndraws);
Avhxs = Avhxs/(imndraws);
Avhxv = Avhxs - Avhxm.^2;
Avhxs = sqrt(Avhxv); % stardard deviation
A0hm = zeros(nvar);
A0hm(a0indx) = Avhxm % mean
A0hv = zeros(nvar);
A0hv(a0indx) = Avhxv; % varaince matrix
A0hs = zeros(nvar);
A0hs(a0indx) = Avhxs; % standar deviation
imfmean = imfmean/(imndraws);
%**** Getting Within-Sequence W and Between-Sequence B_n
% see Gelman, p.331, my Shock(0), pp.12-13, and my Forecast (2), pp.28-31
if (nstarts>1)
AvhxW = (ndraws2/(ndraws2-1))*Avhx_cm/nstarts;
% W: average of j within-sequence variances
AvhxB_n = (nstarts/(nstarts-1)) * ( Avhx_bs/nstarts - (Avhx_bm/nstarts).^2 );
% B/n: variance of J within-sequence means
AvhxB = ndraws2*AvhxB_n; % B
%
B_W1 = AvhxB ./ AvhxW;
B1 = AvhxB;
W1 = AvhxW;
GR1 = sqrt((ndraws2-1)/ndraws2 + B_W1/ndraws2);
% measure of Gelman reduction; need not be 1 to be accurate,
% contrary to what Gelman claims
save outB_W B_W1 B1 W1 GR1 nstarts ndraws2 imndraws timeminutes Avhxs ...
A0xhat A0hm A0hs A0hv IndxGibbs countJump
if ImfErr
save outB_W nswitch imfml imfmean imfcntmulti Range5 ninv invc imstp nvar -append
end
titstr = ['J ' num2str(nstarts) ' n1 ' num2str(ndraws1) ...
' n2 ' num2str(ndraws2) ' timend minutes ' num2str(timend/60)];
disp(' ')
disp(titstr)
disp('B/W sqrt(B) sqrt(W) Std(A0) GR')
format short g
[B_W1 sqrt(B1) sqrt(W1) Avhxs GR1]
else
save outB_W nstarts ndraws1 ndraws2 imndraws timeminutes Avhxs ...
A0xhat A0hm A0hs A0hv IndxGibbs countJump
if ImfErr
save outB_W nswitch imfml imfmean imfcntmulti Range5 ninv invc imstp nvar -append
end
end
|