stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add TZcode submodule to repository

time-shift
Houtan Bastani 2012-04-19 12:03:27 +02:00
parent 9576255ead
commit 6180bc2212
31 changed files with 80 additions and 3050 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitmodules vendored
View File

@ -4,3 +4,6 @@
[submodule "contrib/ms-sbvar/switch_dw"]
path = contrib/ms-sbvar/switch_dw
url = http://www.dynare.org/git/frbatlanta/switch_dw.git
[submodule "contrib/ms-sbvar/TZcode"]
path = contrib/ms-sbvar/TZcode
url = http://www.dynare.org/git/frbatlanta/TZcode.git

1
contrib/ms-sbvar/TZcode Submodule

@ -0,0 +1 @@
Subproject commit 1673c6f7170ac61ad3d3cb832eb3f038c8e4b3ef

80
license.txt
View File

@ -152,13 +152,83 @@ Copyright: 1996 Christopher Sims
2003 Karibzhanov, Waggoner and Zha
License: GPL-3+
Files: matlab/ms-sbvar/cstz/*
Copyright: 1993-2011 Tao Zha
Files: contrib/ms-sbvar/TZcode/*
Copyright: 1997-2012 Tao Zha
License: GPL-3+
Files: matlab/ms-sbvar/cstz/bfgsi.m matlab/ms-sbvar/cstz/csminit.m
matlab/ms-sbvar/cstz/csminwel.m
Copyright: 1993-2011 Tao Zha and Christopher Sims
Files: contrib/ms-sbvar/TZcode/MatlabFiles/szbvar.m
contrib/ms-sbvar/TZcode/MatlabFiles/phg233.m
contrib/ms-sbvar/TZcode/MatlabFiles/phg235.m
contrib/ms-sbvar/TZcode/MatlabFiles/pwf234.m
contrib/ms-sbvar/TZcode/MatlabFiles/pmddf235.m
contrib/ms-sbvar/TZcode/MatlabFiles/pmddf234.m
contrib/ms-sbvar/TZcode/MatlabFiles/pmddf236.m
contrib/ms-sbvar/TZcode/MatlabFiles/csminwel.m
contrib/ms-sbvar/TZcode/MatlabFiles/phg234.m
contrib/ms-sbvar/TZcode/MatlabFiles/pwf235.m
contrib/ms-sbvar/TZcode/MatlabFiles/mnpdf.m
contrib/ms-sbvar/TZcode/MatlabFiles/szasbvar.m
Copyright: 1997-2012 Christopher A. Sims and Tao Zha
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/pmddf233.m
contrib/ms-sbvar/TZcode/MatlabFiles/csminit.m
contrib/ms-sbvar/TZcode/MatlabFiles/gensysoldversion.m
contrib/ms-sbvar/TZcode/MatlabFiles/gensys_z2new.m
contrib/ms-sbvar/TZcode/MatlabFiles/gensysct.m
contrib/ms-sbvar/TZcode/MatlabFiles/csminitworksuntiil0205.m
contrib/ms-sbvar/TZcode/MatlabFiles/gensys_z2.m
contrib/ms-sbvar/TZcode/MatlabFiles/qzdiv.m
contrib/ms-sbvar/TZcode/MatlabFiles/pwf233.m
Copyright: 1997-2012 Christopher A. Sims
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/srestrictrwzalg.m
Copyright: 1997-2012 Juan Rubio-Ramirez, Daniel Waggoner and Tao Zha
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/clmonq.m
contrib/ms-sbvar/TZcode/MatlabFiles/clgls.m
Copyright: 1997-2012 Eric Leeper and Tao Zha
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/a0lhfun.m
Copyright: 1997-2012 Eric Leeper
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/suptitle.m
Copyright: 1997-2012 Drea Thomas
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/fn_gibbsrvar.m
contrib/ms-sbvar/TZcode/MatlabFiles/gibbsvar.m
contrib/ms-sbvar/TZcode/MatlabFiles/simtanzphi.m
contrib/ms-sbvar/TZcode/MatlabFiles/fn_gibbsglb.m
contrib/ms-sbvar/TZcode/MatlabFiles/fn_gibbsrvar_setup.m
contrib/ms-sbvar/TZcode/MatlabFiles/gibbsglb.m
contrib/ms-sbvar/TZcode/MatlabFiles/fn_gibbsrvaroldworks.m
Copyright: 1997-2012 Daniel Waggoner and Tao Zha
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/gensys.m
Copyright: 1996-2012 Christopher A. Sims
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/xydata.m
Copyright: 1997-2012 Lutz Kilian and Tao Zha
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/bfgsi.m
Copyright: 1996 Christopher A. Sims
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/qplot2.m
contrib/ms-sbvar/TZcode/MatlabFiles/ellipse.m
Copyright: 1997-2012 Clark A. Burdick
License: GPL-3+
Files: contrib/ms-sbvar/TZcode/MatlabFiles/chol2.m
Copyright: 1996-2012 Tao Zha
License: GPL-3+
License: GFDL-NIV-1.3+

2
matlab/dynare_config.m
View File

@ -51,8 +51,8 @@ addpath([dynareroot '/kalman/likelihood'])
addpath([dynareroot '/AIM/'])
addpath([dynareroot '/partial_information/'])
addpath([dynareroot '/ms-sbvar/'])
addpath([dynareroot '/ms-sbvar/cstz/'])
addpath([dynareroot '/ms-sbvar/identification/'])
addpath([dynareroot '../contrib/ms-sbvar/TZcode/MatlabFiles/'])
addpath([dynareroot '/parallel/'])
addpath([dynareroot '/particle/'])
addpath([dynareroot '/gsa/'])

46
matlab/ms-sbvar/cstz/bfgsi.m
View File

@ -1,46 +0,0 @@
function H = bfgsi(H0,dg,dx)
% H = bfgsi(H0,dg,dx)
% dg is previous change in gradient; dx is previous change in x;
% 6/8/93 version that updates inverse hessian instead of hessian
% itself.
% Copyright by Christopher Sims 1996. This material may be freely
% reproduced and modified.
dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha)
% Copyright (C) 1996-2011 Tao Zha and Christopher Sims
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if size(dg,2)>1
dg=dg';
end
if size(dx,2)>1
dx=dx';
end
Hdg = H0*dg;
dgdx = dg'*dx;
if (abs(dgdx) >1e-12)
H = H0 + (1+(dg'*Hdg)/dgdx)*(dx*dx')/dgdx - (dx*Hdg'+Hdg*dx')/dgdx;
else
if dispIndx
disp('bfgs update failed.')
disp(['|dg| = ' num2str(sqrt(dg'*dg)) '|dx| = ' num2str(sqrt(dx'*dx))]);
disp(['dg''*dx = ' num2str(dgdx)])
disp(['|H*dg| = ' num2str(Hdg'*Hdg)])
end
H=H0;
end
save H.dat H

216
matlab/ms-sbvar/cstz/csminit.m
View File

@ -1,216 +0,0 @@
function [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,varargin)
% [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,...
% P1,P2,P3,P4,P5,P6,P7,P8)
% retcodes: 0, normal step. 5, largest step still improves too fast.
% 4,2 back and forth adjustment of stepsize didn't finish. 3, smallest
% stepsize still improves too slow. 6, no improvement found. 1, zero
% gradient.
%---------------------
% Modified 7/22/96 to omit variable-length P list, for efficiency and compilation.
% Places where the number of P's need to be altered or the code could be returned to
% its old form are marked with ARGLIST comments.
%
% Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs
% update.
%
% Fixed 7/19/93 to flip eigenvalues of H to get better performance when
% it's not psd.
%
% Fixed 02/19/05 to correct for low angle problems.
%
%tailstr = ')';
%for i=nargin-6:-1:1
% tailstr=[ ',P' num2str(i) tailstr];
%end
% Copyright (C) 1993-2011 Tao Zha and Christopher Sims
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha)
%ANGLE = .03; % when output of this variable becomes negative, we have wrong analytical graident
ANGLE = .005; % works for identified VARs and OLS
%THETA = .03;
THETA = .3; %(0<THETA<.5) THETA near .5 makes long line searches, possibly fewer iterations.
FCHANGE = 1000;
MINLAMB = 1e-9;
% fixed 7/15/94
% MINDX = .0001;
% MINDX = 1e-6;
MINDFAC = .01;
fcount=0;
lambda=1;
xhat=x0;
f=f0;
fhat=f0;
g = g0;
gnorm = norm(g);
%
if (gnorm < 1.e-12) && ~badg % put ~badg 8/4/94
retcode =1;
dxnorm=0;
% gradient convergence
else
% with badg true, we don't try to match rate of improvement to directional
% derivative. We're satisfied just to get some improvement in f.
%
%if(badg)
% dx = -g*FCHANGE/(gnorm*gnorm);
% dxnorm = norm(dx);
% if dxnorm > 1e12
% disp('Bad, small gradient problem.')
% dx = dx*FCHANGE/dxnorm;
% end
%else
% Gauss-Newton step;
%---------- Start of 7/19/93 mod ---------------
%[v d] = eig(H0);
%toc
%d=max(1e-10,abs(diag(d)));
%d=abs(diag(d));
%dx = -(v.*(ones(size(v,1),1)*d'))*(v'*g);
% toc
dx = -H0*g;
% toc
dxnorm = norm(dx);
if dxnorm > 1e12
if dispIndx, disp('Near-singular H problem.'), end
dx = dx*FCHANGE/dxnorm;
end
dfhat = dx'*g0;
%end
%
%
if ~badg
% test for alignment of dx with gradient and fix if necessary
a = -dfhat/(gnorm*dxnorm);
if a<ANGLE
dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g;
% suggested alternate code: ---------------------
dx = dx*dxnorm/norm(dx); % Added 02/19/05 by CAS. This keeps scale invariant to the angle correction
% ------------------------------------------------
dfhat = dx'*g;
% dxnorm = norm(dx); % Removed 02/19/05 by CAS. This line unnecessary with modification that keeps scale invariant
if dispIndx, disp(sprintf('Correct for low angle: %g',a)), end
end
end
if dispIndx, disp(sprintf('Predicted improvement: %18.9f',-dfhat/2)), end
%
% Have OK dx, now adjust length of step (lambda) until min and
% max improvement rate criteria are met.
done=0;
factor=3;
shrink=1;
lambdaMin=0;
lambdaMax=inf;
lambdaPeak=0;
fPeak=f0;
lambdahat=0;
while ~done
if size(x0,2)>1
dxtest=x0+dx'*lambda;
else
dxtest=x0+dx*lambda;
end
% home
f = feval(fcn,dxtest,varargin{:});
%ARGLIST
%f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
% f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8);
if dispIndx, disp(sprintf('lambda = %10.5g; f = %20.7f',lambda,f )), end
%debug
%disp(sprintf('Improvement too great? f0-f: %g, criterion: %g',f0-f,-(1-THETA)*dfhat*lambda))
if f<fhat
fhat=f;
xhat=dxtest;
lambdahat = lambda;
end
fcount=fcount+1;
shrinkSignal = (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | (badg & (f0-f) < 0) ;
growSignal = ~badg & ( (lambda > 0) & (f0-f > -(1-THETA)*dfhat*lambda) );
if shrinkSignal && ( (lambda>lambdaPeak) || (lambda<0) )
if (lambda>0) && ((~shrink) || (lambda/factor <= lambdaPeak))
shrink=1;
factor=factor^.6;
while lambda/factor <= lambdaPeak
factor=factor^.6;
end
%if (abs(lambda)*(factor-1)*dxnorm < MINDX) || (abs(lambda)*(factor-1) < MINLAMB)
if abs(factor-1)<MINDFAC
if abs(lambda)<4
retcode=2;
else
retcode=7;
end
done=1;
end
end
if (lambda<lambdaMax) && (lambda>lambdaPeak)
lambdaMax=lambda;
end
lambda=lambda/factor;
if abs(lambda) < MINLAMB
if (lambda > 0) && (f0 <= fhat)
% try going against gradient, which may be inaccurate
if dispIndx, lambda = -lambda*factor^6, end
else
if lambda < 0
retcode = 6;
else
retcode = 3;
end
done = 1;
end
end
elseif (growSignal && lambda>0) || (shrinkSignal && ((lambda <= lambdaPeak) && (lambda>0)))
if shrink
shrink=0;
factor = factor^.6;
%if ( abs(lambda)*(factor-1)*dxnorm< MINDX ) || ( abs(lambda)*(factor-1)< MINLAMB)
if abs(factor-1)<MINDFAC
if abs(lambda)<4
retcode=4;
else
retcode=7;
end
done=1;
end
end
if ( f<fPeak ) && (lambda>0)
fPeak=f;
lambdaPeak=lambda;
if lambdaMax<=lambdaPeak
lambdaMax=lambdaPeak*factor*factor;
end
end
lambda=lambda*factor;
if abs(lambda) > 1e20;
retcode = 5;
done =1;
end
else
done=1;
if factor < 1.2
retcode=7;
else
retcode=0;
end
end
end
end
if dispIndx, disp(sprintf('Norm of dx %10.5g', dxnorm)), end

306
matlab/ms-sbvar/cstz/csminwel.m
View File

@ -1,306 +0,0 @@
function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,varargin)
%[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit,varargin)
% fcn: string naming the objective function to be minimized
% x0: initial value of the parameter vector
% H0: initial value for the inverse Hessian. Must be positive definite.
% grad: Either a string naming a function that calculates the gradient, or the null matrix.
% If it's null, the program calculates a numerical gradient. In this case fcn must
% be written so that it can take a matrix argument and produce a row vector of values.
% crit: Convergence criterion. Iteration will cease when it proves impossible to improve the
% function value by more than crit.
% nit: Maximum number of iterations.
% varargin: A list of optional length of additional parameters that get handed off to fcn each
% time it is called.
% Note that if the program ends abnormally, it is possible to retrieve the current x,
% f, and H from the files g1.mat and H.mat that are written at each iteration and at each
% hessian update, respectively. (When the routine hits certain kinds of difficulty, it
% write g2.mat and g3.mat as well. If all were written at about the same time, any of them
% may be a decent starting point. One can also start from the one with best function value.)
% NOTE: The display on screen can be turned off by seeting dispIndx=0 in this
% function. This option is used for the loop operation. T. Zha, 2 May 2000
% NOTE: You may want to change stps to 1.0e-02 or 1.0e-03 to get a better convergence. August, 2006
% Copyright (C) 1993-2011 Tao Zha and Christopher Sims
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
Verbose = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha)
dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha)
[nx,no]=size(x0);
nx=max(nx,no);
NumGrad= ( ~ischar(grad) | length(grad)==0);
done=0;
itct=0;
fcount=0;
snit=100;
%tailstr = ')';
%stailstr = [];
% Lines below make the number of Pi's optional. This is inefficient, though, and precludes
% use of the matlab compiler. Without them, we use feval and the number of Pi's must be
% changed with the editor for each application. Places where this is required are marked
% with ARGLIST comments
%for i=nargin-6:-1:1
% tailstr=[ ',P' num2str(i) tailstr];
% stailstr=[' P' num2str(i) stailstr];
%end
f0 = eval([fcn '(x0,varargin{:})']);
%ARGLIST
%f0 = feval(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
% disp('first fcn in csminwel.m ----------------') % Jinill on 9/5/95
if f0 > 1e50, disp('Bad initial parameter.'), return, end
if NumGrad
if length(grad)==0
[g badg] = numgradcd(fcn,x0, varargin{:});
%ARGLIST
%[g badg] = numgradcd(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
else
badg=any(find(grad==0));
g=grad;
end
%numgradcd(fcn,x0,P1,P2,P3,P4);
else
[g badg] = eval([grad '(x0,varargin{:})']);
%ARGLIST
%[g badg] = feval(grad,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
end
retcode3=101;
x=x0;
f=f0;
H=H0;
cliff=0;
while ~done
g1=[]; g2=[]; g3=[];
%addition fj. 7/6/94 for control
if dispIndx
disp('-----------------')
disp('-----------------')
%disp('f and x at the beginning of new iteration')
disp(sprintf('f at the beginning of new iteration, %20.10f',f))
%-----------Comment out this line if the x vector is long----------------
disp([sprintf('x = ') sprintf('%15.8g%15.8g%15.8g%15.8g%15.8g\n',x)]);
end
%-------------------------
itct=itct+1;
[f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,varargin{:});
%ARGLIST
%[f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,P1,P2,P3,P4,P5,P6,P7,...
% P8,P9,P10,P11,P12,P13);
% itct=itct+1;
fcount = fcount+fc;
% erased on 8/4/94
% if (retcode == 1) || (abs(f1-f) < crit)
% done=1;
% end
% if itct > nit
% done = 1;
% retcode = -retcode;
% end
if retcode1 ~= 1
if retcode1==2 || retcode1==4
wall1=1; badg1=1;
else
if NumGrad
[g1 badg1] = numgradcd(fcn, x1,varargin{:});
%ARGLIST
%[g1 badg1] = numgradcd(fcn, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,...
% P10,P11,P12,P13);
else
[g1 badg1] = eval([grad '(x1,varargin{:})']);
%ARGLIST
%[g1 badg1] = feval(grad, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,...
% P10,P11,P12,P13);
end
wall1=badg1;
% g1
save g1.mat g1 x1 f1 varargin;
%ARGLIST
%save g1 g1 x1 f1 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
end
if wall1 % && (~done) by Jinill
% Bad gradient or back and forth on step length. Possibly at
% cliff edge. Try perturbing search direction.
%
%fcliff=fh;xcliff=xh;
if dispIndx
disp(' ')
disp('************************* Random search. *****************************************')
disp('************************* Random search. *****************************************')
disp(' ')
pause(1.0)
end
Hcliff=H+diag(diag(H).*rand(nx,1));
if dispIndx, disp('Cliff. Perturbing search direction.'), end
[f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,varargin{:});
%ARGLIST
%[f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,P1,P2,P3,P4,...
% P5,P6,P7,P8,P9,P10,P11,P12,P13);
fcount = fcount+fc; % put by Jinill
if f2 < f
if retcode2==2 || retcode2==4
wall2=1; badg2=1;
else
if NumGrad
[g2 badg2] = numgradcd(fcn, x2,varargin{:});
%ARGLIST
%[g2 badg2] = numgradcd(fcn, x2,P1,P2,P3,P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
else
[g2 badg2] = eval([grad '(x2,varargin{:})']);
%ARGLIST
%[g2 badg2] = feval(grad,x2,P1,P2,P3,P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
end
wall2=badg2;
% g2
if dispIndx, badg2, end
save g2.mat g2 x2 f2 varargin
%ARGLIST
%save g2 g2 x2 f2 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
end
if wall2
if dispIndx, disp('Cliff again. Try traversing'), end
if norm(x2-x1) < 1e-13
f3=f; x3=x; badg3=1;retcode3=101;
else
gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1);
[f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),varargin{:});
%ARGLIST
%[f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),P1,P2,P3,...
% P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
fcount = fcount+fc; % put by Jinill
if retcode3==2 || retcode3==4
wall3=1; badg3=1;
else
if NumGrad
[g3 badg3] = numgradcd(fcn, x3,varargin{:});
%ARGLIST
%[g3 badg3] = numgradcd(fcn, x3,P1,P2,P3,P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
else
[g3 badg3] = eval([grad '(x3,varargin{:})']);
%ARGLIST
%[g3 badg3] = feval(grad,x3,P1,P2,P3,P4,P5,P6,P7,P8,...
% P9,P10,P11,P12,P13);
end
wall3=badg3;
% g3
if dispIndx, badg3, end
save g3.mat g3 x3 f3 varargin;
%ARGLIST
%save g3 g3 x3 f3 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13;
end
end
else
f3=f; x3=x; badg3=1; retcode3=101;
end
else
f3=f; x3=x; badg3=1;retcode3=101;
end
else
% normal iteration, no walls, or else we're finished here.
f2=f; f3=f; badg2=1; badg3=1; retcode2=101; retcode3=101;
end
else
f1=f; f2=f; f3=f; retcode2=retcode1; retcode3=retcode1;
end
%how to pick gh and xh
if f3<f && badg3==0
if dispIndx, ih=3, end
fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3;
elseif f2<f && badg2==0
if dispIndx, ih=2, end
fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2;
elseif f1<f && badg1==0
if dispIndx, ih=1, end
fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1;
else
[fh,ih] = min([f1,f2,f3]);
if dispIndx, disp(sprintf('ih = %d',ih)), end
%eval(['xh=x' num2str(ih) ';'])
switch ih
case 1
xh=x1;
case 2
xh=x2;
case 3
xh=x3;
end %case
%eval(['gh=g' num2str(ih) ';'])
%eval(['retcodeh=retcode' num2str(ih) ';'])
retcodei=[retcode1,retcode2,retcode3];
retcodeh=retcodei(ih);
if exist('gh')
nogh=isempty(gh);
else
nogh=1;
end
if nogh
if NumGrad
[gh badgh] = feval('numgrad',fcn,xh,varargin{:});
else
[gh badgh] = feval('grad', xh,varargin{:});
end
end
badgh=1;
end
%end of picking
%ih
%fh
%xh
%gh
%badgh
stuck = (abs(fh-f) < crit);
if (~badg)&(~badgh)&(~stuck)
H = bfgsi(H,gh-g,xh-x);
end
if Verbose
if dispIndx
disp('----')
disp(sprintf('Improvement on iteration %d = %18.9f',itct,f-fh))
end
end
if itct > nit
if dispIndx, disp('iteration count termination'), end
done = 1;
elseif stuck
if dispIndx, disp('improvement < crit termination'), end
done = 1;
end
rc=retcodeh;
if rc == 1
if dispIndx, disp('zero gradient'), end
elseif rc == 6
if dispIndx, disp('smallest step still improving too slow, reversed gradient'), end
elseif rc == 5
if dispIndx, disp('largest step still improving too fast'), end
elseif (rc == 4) || (rc==2)
if dispIndx, disp('back and forth on step length never finished'), end
elseif rc == 3
if dispIndx, disp('smallest step still improving too slow'), end
elseif rc == 7
if dispIndx, disp('warning: possible inaccuracy in H matrix'), end
end
f=fh;
x=xh;
g=gh;
badg=badgh;
end
% what about making an m-file of 10 lines including numgrad.m
% since it appears three times in csminwel.m

56
matlab/ms-sbvar/cstz/fn_a0freefun.m
View File

@ -1,56 +0,0 @@
function of = fn_a0freefun(b,Ui,nvar,n0,fss,H0inv)
% of = fn_a0freefun(b,Ui,nvar,n0,fss,H0inv)
%
% Negative logPosterior function for squeesed A0 free parameters, which are b's in the WZ notation
% Note: columns correspond to equations
%
% b: sum(n0)-by-1 vector of A0 free parameters
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation.
% nvar: number of endogeous variables
% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation
% fss: nSample-lags (plus ndobs if dummies are included)
% H0inv: cell(nvar,1). In each cell, posterior inverse of covariance inv(H0) for the ith equation,
% resembling old SpH in the exponent term in posterior of A0, but not divided by T yet.
%----------------
% of: objective function (negative logPosterior)
%
% Tao Zha, February 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
b=b(:); n0=n0(:);
A0 = zeros(nvar);
n0cum = [0;cumsum(n0)];
tra = 0.0;
for kj = 1:nvar
bj = b(n0cum(kj)+1:n0cum(kj+1));
A0(:,kj) = Ui{kj}*bj;
tra = tra + 0.5*bj'*H0inv{kj}*bj; % negative exponential term
end
[A0l,A0u] = lu(A0);
ada = -fss*sum(log(abs(diag(A0u)))); % negative log determinant of A0 raised to power T
of = ada + tra;

59
matlab/ms-sbvar/cstz/fn_a0freegrad.m
View File

@ -1,59 +0,0 @@
function [g,badg] = fn_a0freegrad(b,Ui,nvar,n0,fss,H0inv)
% [g,badg] = a0freegrad(b,Ui,nvar,n0,fss,H0inv)
% Analytical gradient for a0freefun.m in use of csminwel.m. See Dhrymes's book.
%
% b: sum(n0)-by-1 vector of A0 free parameters
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation.
% nvar: number of endogeous variables
% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation
% fss: nSample-lags (plus ndobs if dummies are included)
% H0inv: cell(nvar,1). In each cell, posterior inverse of covariance inv(H0) for the ith equation,
% resembling old SpH in the exponent term in posterior of A0, but not divided by T yet.
%---------------
% g: sum(n0)-by-1 analytical gradient for a0freefun.m
% badg: 0, the value that is used in csminwel.m
%
% Tao Zha, February 2000. Revised, August 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
b=b(:); n0 = n0(:);
A0 = zeros(nvar);
n0cum = [0;cumsum(n0)];
g = zeros(n0cum(end),1);
badg = 0;
%*** The derivative of the exponential term w.r.t. each free paramater
for kj = 1:nvar
bj = b(n0cum(kj)+1:n0cum(kj+1));
g(n0cum(kj)+1:n0cum(kj+1)) = H0inv{kj}*bj;
A0(:,kj) = Ui{kj}*bj;
end
B=inv(A0');
%*** Add the derivative of -Tlog|A0| w.r.t. each free paramater
for ki = 1:sum(n0)
n = max(find( (ki-n0cum)>0 )); % note, 1<=n<=nvar
g(ki) = g(ki) - fss*B(:,n)'*Ui{n}(:,ki-n0cum(n));
end

75
matlab/ms-sbvar/cstz/fn_calyrqm.m
View File

@ -1,75 +0,0 @@
function [Myrqm,nMyrqm] = fn_calyrqm(q_m,Byrqm,Eyrqm)
% [Myrqm,nMyrqm] = fn_calyrqm(q_m,Byrqm,Eyrqm)
%
% Given the beginning and end years and quarters (months), export a matrix of all years and
% quarters (months) for these years and in between
%
% q_m: 4 if quarterly and 12 if monthly
% Byrqm: [year quarter(month)] -- all integers, the begining year and quarter (month)
% Eyrqm: [year quarter(month)] -- all integers, the end year and quarter (month)
%-------------------
% Myrqm: matrix of all years and quarters (months) between and incl. Byrqm and Eyrqm
% nMyrqm: number of data points incl. Byrqm and Eyrqm
%
% Tao Zha, April 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if ~isempty(find(Byrqm-round(Byrqm))) | (q_m-round(q_m)) | ~isempty(find(Byrqm-round(Byrqm)))
error('argin qm, Byrqm, or Eyrqm must of integer')
elseif Byrqm(1)>Eyrqm(1)
error('Eyrqm(1) must be equal to or greater than Byrqm(1)')
elseif Byrqm(1)==Eyrqm(1)
if Byrqm(2)>Eyrqm(2)
error('Eyrqm(2) must be equal to or greater than Byrqm(2) because of the same year')
end
end
Yr = Byrqm(1)+[0:Eyrqm(1)-Byrqm(1)]';
if length(Yr)>=2 % there are years and quarters (months) between Byrqm and Eyrqm
n=length(Yr)-2;
C=zeros(n*q_m,2);
C(:,1) = kron(Yr(2:end-1),ones(q_m,1));
C(:,2) = kron(ones(n,1),[1:q_m]');
%* initialize a matrix of years and quarters (months) including Byrqm and Eyrqm
Myrqm = zeros((q_m-Byrqm(2)+1)+Eyrqm(2)+n*q_m,2);
%* Years in between
n1=q_m-Byrqm(2)+1;
n2=Eyrqm(2);
Myrqm(n1+1:end-n2,:) = C;
%* Beginning year
for k=1:n1
Myrqm(k,:) = [Byrqm(1) Byrqm(2)+k-1];
end
%* End year
for k=1:n2
Myrqm(end-Eyrqm(2)+k,:) = [Eyrqm(1) k];
end
else %* all the data are in the same calendar year
n1=Eyrqm(2)-Byrqm(2)+1;
Myrqm = zeros(n1,2);
for k=1:n1
Myrqm(k,:) = [Byrqm(1) Byrqm(2)+k-1];
end
end
nMyrqm = size(Myrqm,1);

60
matlab/ms-sbvar/cstz/fn_dataext.m
View File

@ -1,60 +0,0 @@
function [xdsube,Brow,Erow] = fn_dataext(Byrqm,Eyrqm,xdatae)
% xdsube = dataext(Byrqm,Eyrqm,xdatae)
% Extract subset of xdatae, given Byrqm and Eyrqm
%
% Byrqm: [year quarter(month)]: beginning year and period. If Byqm(2)=0, we get annual data.
% Eyrqm: [year period]: end year and period. If Byrqm(2)=0, it must be Eyrqm(2)=0.
% xdatae: all data (some of which may be NaN) with the first 2 columns indicating years and periods.
%----------
% xdsube: subset of xdatae.
% Brow: the row number in xdatee that marks the first row of xdsube.
% Erow: the row number in xdatee that marks the last row of xdsube.
%
% Tao Zha, April 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if (Byrqm(2)==0) && (Eyrqm(2)~=0)
error('If annual data, make sure both Byrqm(2) and Eyrqm(2) are zero')
end
Brow = min(find(xdatae(:,1)==Byrqm(1)));
if isempty(Brow)
error('Byrqm is outside the date range indicated by xdatae(:,1:2)!')
end
if Byrqm(2)>0
nadt=Byrqm(2)-xdatae(Brow,2);
if nadt<0
error('Byrqm is outside the date range indicated by xdatae(:,1:2)!')
end
Brow=Brow+nadt;
end
%
Erow = min(find(xdatae(:,1)==Eyrqm(1)));
if isempty(Erow)
error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!')
end
if Eyrqm(2)>0
nadt2=Eyrqm(2)-xdatae(Erow,2);
if nadt<0
error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!')
end
Erow=Erow+nadt2;
end
%
xdsube = xdatae(Brow:Erow,:); % with dates

117
matlab/ms-sbvar/cstz/fn_datana.m
View File

@ -1,117 +0,0 @@
function [yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(xdatae,q_m,vlistlog,vlistper,Byrqm,Eyrqm)
% [yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(Byrqm,Eyrqm,xdatae,q_m,vlistlog,vlistper)
%
% Generate prior period, period-to-period-in-last-year, and annual growth rates
% For annual rates, works for both calendar and any annual years, depending on Byrqm and Eyrqm
% If monthly data, we haven't got time to get quarterly growth rates yet.
%
% xdatae: all data (logged levels or interest rates/100, some of which may be NaN) with the first
% 2 columns indicating years and periods.
% q_m: quarter or month period
% vlistlog: sublist for logged variables
% vlistper: sublists for percent variables
% Byrqm: [year quarter(month)]: beginning year and period. If Byqm(2)~=1, we don't get
% calendar annual rates. In other words, the first column of yactyge (which
% indicates years) does not mean calendar years. Byqm(2) must be specified; in other
% words, it must be not set to 0 as in, say, fn_dataext.
% Eyrqm: [year period]: end year and period. Eyqm(2) must be specified; in other words, it
% must be not set to 0 as in, say, fn_dataext.
% NOTE: if no inputs Byrqm and Eyrqm are specified, all growth rates begin at xdatae(1,1:2).
%----------
% yactyrge: annual growth rates with dates in the first 2 columns.
% yactyre: annual average logged level with dates in the 1st 2 columns.
% yactqmyge: period-to-period-in-last-year annual growth rates with dates in the first 2 columns.
% yactqmge: prior-period annualized growth rates with dates in the first 2 columns.
% yactqme: data (logged levels or interest rates/100) with dates in the first 2 columns.
% Same as xdatae but with Brow:Erow.
%
% Tao Zha, April 2000.
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if size(xdatae,1)<2*q_m
error('We need at least two years of xdatae to get annual rates. Check xdatae!!')
end
if nargin==4
Brow=1; Erow=length(xdatae(:,1));
nyr = floor((Erow-Brow+1)/q_m);
yrsg = [xdatae(q_m+1,1):xdatae(q_m+1,1)+nyr-2]'; % for annual growth later on
else
if Byrqm(2)<1 || Eyrqm(2)<1
error('This function requires specifying both years and months (quarters) in Byrqm and Eyrqm')
end
Brow = min(find(xdatae(:,1)==Byrqm(1)));
if isempty(Brow)
error('Byrqm is outside the date range of xdatae(:,1:2)!')
end
nadt=Byrqm(2)-xdatae(Brow,2);
if nadt<0
error('Byrqm is outside the date range indicated by xdatae(:,1:2)!')
end
Brow=Brow+nadt;
%
Erow = min(find(xdatae(:,1)==Eyrqm(1)));
if isempty(Brow)
error('Eyrqm is outside the date range of xdatae(:,1:2)!')
end
nadt2=Eyrqm(2)-xdatae(Erow,2);
if nadt<0
error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!')
end
Erow=Erow+nadt2;
nyr = floor((Erow-Brow+1)/q_m);
yrsg = [Byrqm(1)+1:Byrqm(1)+nyr-1]'; % for annual growth later on, which will
% start at Byrqm(1) instead of Byrqm(1)+1
end
%
yactqme = xdatae(Brow:Erow,:); % with dates
yactqm = yactqme(:,3:end); % only data
%======== prior period change (annaluized rate)
yactqmg = yactqm(2:end,:); % start at second period to get growth rate
yactqmg(:,vlistlog) = (yactqm(2:end,vlistlog) - yactqm(1:end-1,vlistlog)) .* q_m;
% monthly, 12*log(1+growth rate), annualized growth rate
yactqmg(:,vlistlog) = 100*(exp(yactqmg(:,vlistlog))-1);
yactqmg(:,vlistper) = 100*yactqmg(:,vlistper);
yactqmge = [yactqme(2:end,1:2) yactqmg];
%======== change from the last year
yactqmyg = yactqm(q_m+1:end,:); % start at the last-year period to get growth rate
yactqmyg(:,vlistlog) = (yactqm(q_m+1:end,vlistlog) - yactqm(1:end-q_m,vlistlog));
yactqmyg(:,vlistlog) = 100*(exp(yactqmyg(:,vlistlog))-1);
yactqmyg(:,vlistper) = 100*yactqmyg(:,vlistper);
yactqmyge = [yactqme(q_m+1:end,1:2) yactqmyg];
%======== annual growth rates
nvar = length(xdatae(1,3:end));
ygmts = yactqm(1:nyr*q_m,:); % converted to the multiplication of q_m
ygmts1 = reshape(ygmts,q_m,nyr,nvar);
ygmts2 = sum(ygmts1,1) ./ q_m;
ygmts3 = reshape(ygmts2,nyr,nvar); % converted to annual average series
%
yactyrg = ygmts3(2:end,:); % start at the last-year period to get growth rate
yactyrg(:,vlistlog) = ygmts3(2:end,vlistlog) - ygmts3(1:end-1,vlistlog);
% annual rate: log(1+growth rate)
yactyrg(:,vlistlog) = 100*(exp(yactyrg(:,vlistlog))-1);
yactyrg(:,vlistper) = 100*yactyrg(:,vlistper);
yactyrge = [yrsg zeros(nyr-1,1) yactyrg];
yrsg1=[yrsg(1)-1:yrsg(end)]';
yactyre = [yrsg1 zeros(nyr,1) ygmts3];

164
matlab/ms-sbvar/cstz/fn_dataxy.m
View File

@ -1,164 +0,0 @@
function [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh,e] = fn_dataxy(nvar,lags,z,mu,indxDummy,nexo)
% [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh,e] = fn_dataxy(nvar,lags,z,mu,indxDummy,nexo)
%
% Export arranged data matrices for future estimation of the VAR.
% If indxDummy=0, no mu's are used at all and Bh is the OLS estimate.
% If indxDummy=1, only mu(5) and mu(6) as dummy observations. See fn_rnrprior*.m for using mu(1)-mu(4).
% See Wagonner and Zha's Gibbs sampling paper.
%
% nvar: number of endogenous variables.
% lags: the maximum length of lag
% z: T*(nvar+(nexo-1)) matrix of raw or original data (no manipulation involved)
% with sample size including lags and with exogenous variables other than a constant.
% Order of columns: (1) nvar endogenous variables; (2) (nexo-1) exogenous variables;
% (3) constants will be automatically put in the last column.
% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where
% only mu(5) and mu(6) are used for dummy observations in this function (i.e.,
% mu(1)-mu(4) are irrelevant here). See fn_rnrprior*.m for using mu(1)-mu(4).
% mu(1): overall tightness and also for A0; (0.57)
% mu(2): relative tightness for A+; (0.13)
% mu(3): relative tightness for the constant term; (0.1)
% mu(4): tightness on lag decay; (1)
% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5)
% mu(6): weight on single dummy initial observation including constant
% (cointegration, unit roots, and stationarity); (5)
% indxDummy: 1: add dummy observations to the data; 0: no dummy added.
% nexo: number of exogenous variables. The constant term is the default setting. Besides this term,
% we have nexo-1 exogenous variables. Optional. If left blank, nexo is set to 1.
% -------------------
% xtx: X'X: k-by-k where k=ncoef
% xty: X'Y: k-by-nvar
% yty: Y'Y: nvar-by-nvar
% fss: T: sample size excluding lags. With dummyies, fss=nSample-lags+ndobs.
% phi: X; T-by-k; column: [nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term]
% y: Y: T-by-nvar where T=fss
% ncoef: number of coefficients in *each* equation. RHS coefficients only, nvar*lags+nexo
% xr: the economy size (ncoef-by-ncoef) in qr(phi) so that xr=chol(X'*X) or xr'*xr=X'*X
% Bh: ncoef-by-nvar estimated reduced-form parameter; column: nvar;
% row: ncoef=[nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term]
% e: estimated residual e = y -x*Bh, T-by-nvar
%
% Tao Zha, February 2000.
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin == 5
nexo=1; % default for constant term
elseif nexo<1
error('We need at least one exogenous term so nexo must >= 1')
end
%*** original sample dimension without dummy prior
nSample = size(z,1); % the sample size (including lags, of course)
sb = lags+1; % original beginning without dummies
ncoef = nvar*lags+nexo; % number of coefficients in *each* equation, RHS coefficients only.
if indxDummy % prior dummy prior
%*** expanded sample dimension by dummy prior
ndobs=nvar+1; % number of dummy observations
fss = nSample+ndobs-lags;
%
% **** nvar prior dummy observations with the sum of coefficients
% ** construct X for Y = X*B + U where phi = X: (T-lags)*k, Y: (T-lags)*nvar
% ** columns: k = # of [nvar for 1st lag, ..., nvar for last lag, exo var, const]
% ** Now, T=T+ndobs -- added with "ndobs" dummy observations
%
phi = zeros(fss,ncoef);
%* constant term
const = ones(fss,1);
const(1:nvar) = zeros(nvar,1);
phi(:,ncoef) = const; % the first nvar periods: no or zero constant!
%* other exogenous (than) constant term
phi(ndobs+1:end,ncoef-nexo+1:ncoef-1) = z(lags+1:end,nvar+1:nvar+nexo-1);
exox = zeros(ndobs,nexo);
phi(1:ndobs,ncoef-nexo+1:ncoef-1) = exox(:,1:nexo-1);
% this = [] when nexo=1 (no other exogenous than constant)
xdgel = z(:,1:nvar); % endogenous variable matrix
xdgelint = mean(xdgel(1:lags,:),1); % mean of the first lags initial conditions
%* Dummies
for k=1:nvar
for m=1:lags
phi(ndobs,nvar*(m-1)+k) = xdgelint(k);
phi(k,nvar*(m-1)+k) = xdgelint(k);
% <<>> multiply hyperparameter later
end
end
%* True data
for k=1:lags
phi(ndobs+1:fss,nvar*(k-1)+1:nvar*k) = xdgel(sb-k:nSample-k,:);
% row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, exo var, const]
% Thus, # of columns is nvar*lags+nexo = ncoef.
end
%
% ** Y with "ndobs" dummies added
y = zeros(fss,nvar);
%* Dummies
for k=1:nvar
y(ndobs,k) = xdgelint(k);
y(k,k) = xdgelint(k);
% multiply hyperparameter later
end
%* True data
y(ndobs+1:fss,:) = xdgel(sb:nSample,:);
phi(1:nvar,:) = 1*mu(5)*phi(1:nvar,:); % standard Sims and Zha prior
y(1:nvar,:) = mu(5)*y(1:nvar,:); % standard Sims and Zha prior
phi(nvar+1,:) = mu(6)*phi(nvar+1,:);
y(nvar+1,:) = mu(6)*y(nvar+1,:);
[xq,xr]=qr(phi,0);
xtx=xr'*xr;
xty=phi'*y;
[yq,yr]=qr(y,0);
yty=yr'*yr;
Bh = xr\(xr'\xty); % xtx\xty where inv(X'X)*(X'Y)
e=y-phi*Bh;
else
fss = nSample-lags;
%
% ** construct X for Y = X*B + U where phi = X: (T-lags)*k, Y: (T-lags)*nvar
% ** columns: k = # of [nvar for 1st lag, ..., nvar for last lag, exo var, const]
%
phi = zeros(fss,ncoef);
%* constant term
const = ones(fss,1);
phi(:,ncoef) = const; % the first nvar periods: no or zero constant!
%* other exogenous (than) constant term
phi(:,ncoef-nexo+1:ncoef-1) = z(lags+1:end,nvar+1:nvar+nexo-1);
% this = [] when nexo=1 (no other exogenous than constant)
xdgel = z(:,1:nvar); % endogenous variable matrix
%* True data
for k=1:lags
phi(:,nvar*(k-1)+1:nvar*k) = xdgel(sb-k:nSample-k,:);
% row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, exo var, const]
% Thus, # of columns is nvar*lags+nexo = ncoef.
end
%
y = xdgel(sb:nSample,:);
[xq,xr]=qr(phi,0);
xtx=xr'*xr;
xty=phi'*y;
[yq,yr]=qr(y,0);
yty=yr'*yr;
Bh = xr\(xr'\xty); % xtx\xty where inv(X'X)*(X'Y)
e=y-phi*Bh;
end

33
matlab/ms-sbvar/cstz/fn_ergodp.m
View File

@ -1,33 +0,0 @@
function gpi = fn_ergodp(P)
% gpi = fn_ergodp(P)
% Compute the ergodic probabilities. See Hamilton p.681.
%
% P: n-by-n matrix of transition matrix where all elements in each column sum up to 1.
%-----
% gpi: n-by-1 vector of ergodic probabilities.
%
% Tao Zha August 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
[gpim,gpid] = eig(P); % m: matrix; d: diagonal
[gpidv,gpidvinx] = sort(diag(gpid));
gpidv = fliplr(gpidv);
gpidvinx = flipud(gpidvinx);
gpim = gpim(:,gpidvinx);
gpi = gpim(:,1)/sum(gpim(:,1));

322
matlab/ms-sbvar/cstz/fn_fcstidcnd.m
View File

@ -1,322 +0,0 @@
function [yhat,Estr,rcon,Rcon,u,v,d] = fn_fcstidcnd(valuecon,stepcon,varcon,nstepsm,...
nconstr,eq_ms,nvar,lags,phil,Aband,Sband,yfore_h,imf3s_h,A0_h,Bh_h,...
forep,TLindx,TLnumber,nCms,eq_Cms)
% [yhat,Estr,rcon,Rcon,u,v,d] = fn_fcstidcnd(valuecon,stepcon,varcon,nstepsm,...
% nconstr,eq_ms,nvar,lags,phil,Aband,Sband,yfore_h,imf3s_h,A0_h,Bh_h,...
% forep,TLindx,TLnumber,nCms,eq_Cms)
%
% Conditional forecasting in the identified model with or without error bands
% It handles conditions on average values as well, so "valuecon" must be
% expressed at average (NOT sum) level.
% Aband is used only once when nconstr>0 and Aband=1, where Gibbs sampler may be used
% Unconditional forecast when imf3s_h, etc is fixed and nconstr=0.
%
% valuecon: vector of values conditioned
% stepcon: sequence (cell) of steps conditioned; if length(stepcon{i}) > 1, the condition
% is then an arithmetic average of log(y) over the stepcon{i} period.
% varcon: vector of variables conditioned
% nconstr: number of DLS constraints
% nstepsm: maximum number of steps in all DLS constraints
% nvar: number of variables in the BVAR model
% lags: number of lags in the BVAR model
% phil: the 1-by-(nvar*lags+1) data matrix where k=nvar*lags+1
% (last period plus lags before the beginning of forecast)
% Aband: 1: draws from A0 and Bh; 0: no draws
% Sband: 1: draws from random shocks E; 0: no random shocks
% yfore_h: uncondtional forecasts: forep-by-nvar. Never used when nconstr=0.
% In this case, set it to [];
% imf3s_h: 3-dimensional impulse responses matrix: impsteps-by-nvar shocks-by-nvar responses
% Never used when nconstr=0. In this case, set it to [];
% A0_h: A0 contemporaneous parameter matrix
% Bh_h: reduced-form parameter matrix: k-by-nvar, y(t) = X(t)*Bh+e(t)
% where X(t) is k-by-nvar and y(t) is 1-by-nvar
% forep: # of forecast periods (e.g., monthly for a monthly model)
% TLindx: 1-by-nCms vector of 1's and 0's, indicating tight or loose; 1: tighter, 0: looser
% Used only when /* (MS draws) is activated. Right now, MS shocks are deterministic.
% TLnumber: 1-by-nCms vector, lower bound for tight and upper bound for loose
% nCms: # of LZ conditions
% eq_Cms: equation location of MS shocks
% ------
% yhat: conditional forecasts: forep-by-nvar
% Estr: backed-out structural shocks (from N(0,1))
% rcon: vector - the difference between valuecon and log(yfore) (unconditional forecasts)
% Rcon: k-by-q (q constranits and k=nvar*max(nsteps)) so that
% Rcon'*e = rcon where e is k-by-1
% [u,d,v]: svd(Rcon,0)
%
%% See Zha's note "Forecast (1)" p. 5, RATS manual (some errors in RATS), etc.
%
%% Some notations: y(t+1) = y(t)B1 + e(t+1)inv(A0). e(t+1) is 1-by-n.
%% Let r(t+1)=e(t+1)inv(A0) + e(t+2)C + .... where inv(A0) is impulse
%% response at t=1, C at t=2, etc. The row of inv(A0) or C is
%% all responses to one shock.
%% Let r be q-by-1 (such as r(1) = r(t+1)
%% = y(t+1) (constrained) - y(t+1) (forecast)).
%% Use impulse responses to find out R (k-by-q) where k=nvar*nsteps
%% where nsteps the largest constrained step. The key of the program
%% is to creat R using impulse responses
%% Optimal solution for shock e where R'*e=r and e is k-by-1 is
%% e = R*inv(R'*R)*r and k>=q
%
% Copyright (c) March 1998 by Tao Zha. Revised November 1998;
% 3/20/99 Disenabled draws of MS shcoks. To enable it, activate /* part
% 3/20/99 Added A0_h and forep and deleted Cms as input argument. Previous
% programs may not be compatible.
% 3/15/2004 There are some BUG problems when calling fn_fcstcnd.m().
% Copyright (C) 1998-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
DLSIdShock = ~isempty(eq_ms); % if not empty, the MS shock is identified as in DLS
impsteps=size(imf3s_h,1);
if (forep<nstepsm) || (impsteps<nstepsm)
disp('Increase # of forecast or impulse steps!!')
disp('Or decrease # of constraints (nconstr) or constrained steps (stepcon(i))!!')
error('Maximum of conditional steps > # of forecast or impulse steps!!')
end
kts = nvar*nstepsm; % k -- ts: total shocks some of which are restricted and others
% are free.
%*** initializing
Rcon = zeros(kts,nconstr); % R: k-by-q
Econ = zeros(kts,1); % E: k-by-1
rcon = zeros(nconstr,1); % r: q-by-1
%rcon=valuecon-diag(yfore(stepcon,varcon)); % another way is to use "loop" below.
tcwc = nvar*lags; % total coefficients without constant
phi=phil;
%----------------------------------------------------
% Form rcon, Rcon, and Econ (the mean of structural shocks)
%----------------------------------------------------
if nconstr
A0in = reshape(imf3s_h(1,:,:),nvar,nvar); % nvar shocks-by-nvar responses
for i=1:nconstr
rcon(i)=length(stepcon{i})*valuecon(i) - ...
sum(yfore_h(stepcon{i},varcon(i)),1); %<<>>
Rmat = zeros(nstepsm,nvar);
r2mat = zeros(nstepsm,1); % simply one identified equation
% Must be here inside the loop because it's matrix of one column of Rcon
for j=1:length(stepcon{i})
if DLSIdShock % Assuming the Fed can't see all other shocks within a month
Rmat(1:stepcon{i}(j),eq_ms) = Rmat(1:stepcon{i}(j),eq_ms) + ...
imf3s_h(stepcon{i}(j):-1:1,eq_ms,varcon(i));
% Rmat: row--nstepsm, column--nvar shocks (here all shocks except
% the identified one are set to zero) for a particular
% endogenous variable 'varcon(i)'. See Zha Forcast (1), pp.6-7
else % Rcon random with (A0,A+)
Rmat(1:stepcon{i}(j),:) = Rmat(1:stepcon{i}(j),:) + ...
imf3s_h(stepcon{i}(j):-1:1,:,varcon(i));
% Rmat: row--nstepsm, column--nvar shocks (here all shocks are
% *not* set to zero) for a particular endogenous
% variable 'varcon(i)'. See Zha Forcast (1), pp.6-7
end
end
Rmatt = Rmat'; % Now, nvar-by-nstepsm. I think here is where RATS has an error
% i.e. "OVERR" is not transposed when overlaid to "CAPR"
Rcon(:,i)=Rmatt(:); % Rcon: k-by-q where q=nconstr
end
[u d v]=svd(Rcon,0); %trial
%???? Can we reduce the time by computing inv(R'*R) directly?
% rtr = Rcon'*Rcon; %trial
% rtrinv = inv(Rcon'*Rcon); %trial
vd=v.*(ones(size(v,2),1)*diag(d)'); %trial
dinv = 1./diag(d); % inv(diag(d))
vdinv=v.*(ones(size(v,2),1)*dinv'); %trial
rtr=vd*vd'; % R'*R
rtrinv = vdinv*vdinv'; % inv(R'*R)
Econ = Rcon*rtrinv*rcon; % E = R*inv(R'R)*r; the mean of structural shocks
else
Econ = zeros(kts,1); % the mean of shocks is zero under no variable condition
Rcon = NaN;
rcon = NaN;
u = NaN;
d = NaN;
v = NaN;
end
%---------------------------------------
% No uncertainty at all or only random (A0,A+)
% In other words, no future shocks
%---------------------------------------
if (~Sband) %|| (nconstr && (length(eq_ms)==1))
% length(eq_ms)==1 implies one-one mapping between MS shocks and, say, FFR
% if nstepsm==nconstr. If this condition does not hold, this procedure
% is incorrect. I don't have time to fix it now (3/20/99). So I use
% this as a proximation
Estr = reshape(Econ,nvar,nstepsm);
Estr = Estr'; % transpose so that
% Estr: structural shocks. Row--steps, Column--n shocks
Estr = [Estr;zeros(forep-nstepsm,nvar)];
% Now, forep-by-nvar -- ready for forecasts
Estr(1:nCms,eq_Cms) = TLnumber(:);
Ures = Estr/A0_h; % nstepsm-by-nvar
% Ures: reduced-form residuals. Row--steps; Column--n shocks
% ** reconstruct x(t) for y(t+h) = x(t+h-1)*B
% ** where phi = x(t+h-1) with last column being constant
%
yhat = zeros(forep,nvar);
for k=1:forep
yhat(k,:) = phi*Bh_h + Ures(k,:);
phi(1,nvar+1:tcwc) = phi(1,1:tcwc-nvar);
phi(1,1:nvar) = yhat(k,:);
%
end
%---------------------------------------
% With random future shocks and possibly (A0,A+) depending
% on if imf3s_h is random
%---------------------------------------
else
%--------------
% Condition on variables and A random
%--------------
if nconstr && Aband
warning(' ')
disp('This situation (both E and A random) is still under construction')
disp('It is closely related to Waggoner and Zha ReStat Gibbs sampling method')
disp('Please press ctrl-c to abort')
pause
elseif nconstr
%--------------
% Condition on variables and DLS MS shock, no A random but S random
%--------------
if DLSIdShock % other shocks are indepedent of the eq_ms shock
% 3/20/99 The following may be problematic because Osk should depend
% on u (A0_h and Bh_h) in general. I haven't worked out any good version
%/*
% Osk = randn(kts,1); % other shocks
% for j=1:nstepsm
% Osk(nvar*(j-1)+eq_ms)=0; % no shock to the MS or identified equation
% end
% Estr = Econ + Osk; % Econ is non zero only at position
% % eq_ms*j where j=1:nstepsm
% Estr = reshape(Estr,nvar,nstepsm);
% Estr = Estr'; % transpose so that
% % Estr: structural shocks. Row--steps, Column--n shocks
% Estr = [Estr;randn(forep-nstepsm,nvar)];
% % Now, forep-by-nvar -- ready for forecasts
%
disp('DLS')
Ome = eye(kts) - u*u'; % note, I-u*u' = I - R*inv(R'*R)*R'
%[u1 d1 v1] = svd(Ome); % too slow
[u1 d1] = eig(Ome);
Stdcon = u1*diag(sqrt(diag(abs(d1)))); % lower triagular chol of conditional variance
% see Zha's forecast (1), p.17
tmp1=zeros(nvar,nstepsm);
tmp1(eq_ms,:)=randn(1,nstepsm);
tmp2=tmp1(:);
%Estr1 = Econ + Stdcon*randn(kts,1);
%jnk = reshape(Stdcon*tmp2,nvar,nstepsm)
Estr1 = Econ + Stdcon*tmp2;
Estr2 = reshape(Estr1,nvar,nstepsm);
Estr2 = Estr2'; % transpose so that
% Estr2: structural shocks. Row--nstepsm, Column--n shocks
Estr = [Estr2;randn(forep-nstepsm,nvar)];
% Now, forep-by-nvar -- ready for forecasts
else
Ome = eye(kts) - u*u'; % note, I-u*u' = I - R*inv(R'*R)*R'
%[u1 d1 v1] = svd(Ome); % too slow
[u1 d1] = eig(Ome);
Stdcon = u1*diag(sqrt(diag(abs(d1)))); % lower triagular chol of conditional variance
% see Zha's forecast (1), p.17
%--------------
% Condition on variables and LZ MS shock, no A random but S random
% This section has not be tested yet, 10/14/98
%--------------
if nCms
Estr1 = Econ + Stdcon*randn(kts,1);
Estr2 = reshape(Estr1,nvar,nstepsm);
Estr2 = Estr2'; % transpose so that
% Estr2: structural shocks. Row--nstepsm, Column--n shocks
Estr = [Estr2;randn(forep-nstepsm,nvar)];
% Now, forep-by-nvar -- ready for forecasts
Estr(1:nCms,eq_Cms) = TLnumber(:);
%/* draw MS shocks
% for k=1:nCms
% if TLindx(k) % tighter
% while (Estr(k,eq_Cms)<TLnumber(k))
% Estr(k,eq_Cms) = randn(1,1);
% end
% else % looser
% while (Estr(k,eq_Cms)>TLnumber(k))
% Estr(k,eq_Cms) = randn(1,1);
% end
% end
% end
%--------------
% Condition on variables only, no A random but S random
%--------------
else
Estr1 = Econ + Stdcon*randn(kts,1);
Estr2 = reshape(Estr1,nvar,nstepsm);
Estr2 = Estr2'; % transpose so that
% Estr2: structural shocks. Row--nstepsm, Column--n shocks
Estr = [Estr2;randn(forep-nstepsm,nvar)];
% Now, forep-by-nvar -- ready for forecasts
end
end
%--------------
% Condition on LZ MS shocks only, S random and possibly A random depending on
% if A0_h and Bh_h are random
%--------------
else
if nCms
Estr = randn(forep,nvar);
% Now, forep-by-nvar -- ready for forecasts
Estr(1:nCms,eq_Cms) = TLnumber(:);
%/* draw MS shocks
% for k=1:nCms
% if TLindx(k) % tighter
% while (Estr(k,eq_Cms)<TLnumber(k))
% Estr(k,eq_Cms) = randn(1,1);
% end
% else % looser
% while (Estr(k,eq_Cms)>TLnumber(k))
% Estr(k,eq_Cms) = randn(1,1);
% end
% end
% end
else
Estr = randn(forep,nvar); % Unconditional forecast
% Now, forep-by-nvar -- ready for forecasts
end
end
%
Ures = Estr/A0_h; % nstepsm-by-nvar
% Ures: reduced-form residuals. Row--steps; Column--n shocks
% ** reconstruct x(t) for y(t+h) = x(t+h-1)*B
% ** where phi = x(t+h-1) with last column being constant
%
yhat = zeros(forep,nvar);
for k=1:forep
yhat(k,:) = phi*Bh_h + Ures(k,:);
phi(1,nvar+1:tcwc) = phi(1,1:tcwc-nvar);
phi(1,1:nvar) = yhat(k,:);
end
end

74
matlab/ms-sbvar/cstz/fn_forecast.m
View File

@ -1,74 +0,0 @@
function yhat = fn_forecast(Bh,phi,nn,nexo,Xfexo)
% yhat = fn_forecast(Bh,phi,nn,nexo,Xfexo)
% Unconditional forecating without shocks.
% y_hat(t+h) = c + x_hat(t+h-1)*Bh, X: 1*k; Bh: k*nvar; y_hat: 1*nvar
%
% Bh: k-by-nvar, the (posterior) estimate of B.
% phi: the 1-by-(nvar*lags+nexo) data matrix X where k=nvar*lags+1
% (last period plus lags before the beginning of forecast).
% nn: [nvar,lags,nfqm], nfqm: forecast periods (months or quarters).
% nexo: number of exogenous variables. The constant term is the default setting.
% Besides this term, we have nexo-1 exogenous variables.
% Xfexo: nfqm-by-nexo-1 vector of exoenous variables in the forecast horizon where
% nfqm: number of forecast periods.
%-----------
% yhat: nfqm-by-nvar forecast.
%
% See fn_forecastsim.m with shocks; fn_forecaststre.m.
% Copyright (C) 2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin == 3
nexo=1; % default for constant term
elseif nexo<1
error('We need at least one exogenous term so nexo must >= 1')
end
% ** setup
nvar = nn(1);
lags = nn(2);
nfqm = nn(3);
tcwx = nvar*lags; % total coefficeint without exogenous variables
if nexo>1
if (nfqm > size(Xfexo,1))
disp(' ')
warning('Make sure the forecast horizon in the exogenous variable matrix Xfexo > forecast periods')
disp('Press ctrl-c to abort')
pause
elseif ((nexo-1) ~= size(Xfexo,2))
disp(' ')
warning('Make sure that nexo matchs the exogenous variable matrix Xfexo')
disp('Press ctrl-c to abort')
pause
end
end
% ** reconstruct x(t) for y(t+h) = x(t+h-1)*B
% ** where phi = x(t+h-1) with last column being constant
yhat = zeros(nfqm,nvar);
for k=1:nfqm
yhat(k,:) = phi*Bh;
%*** lagged endogenous variables
phi(1,nvar+1:tcwx) = phi(1,1:tcwx-nvar);
phi(1,1:nvar) = yhat(k,:);
%*** exogenous variables excluding constant terms
if (nexo>1)
phi(1,tcwx+1:tcwx+nexo-1) = Xfexo(k,:);
end
end

65
matlab/ms-sbvar/cstz/fn_foregraph.m
View File

@ -1,65 +0,0 @@
function fn_foregraph(yfore,yacte,keyindx,rnum,cnum,q_m,ylab,forelabel,conlab)
%
% Graph annual (or calendar year) forecasts vs actual (all data from "msstart.m")
%
% yfore: actual and forecast annual growth data with dates.
% yacte: actual annual growth data with dates.
% keyindx: index for the variables to be graphed
% rnum: number of rows in subplot
% cnum: number of columns in subplot
% q_m: if 4 or 12, quarterly or monthly data
% ylab: string array for the length(keyindx)-by-1 variables
% forelabel: title label for as of time of forecast
% conlab: label for what conditions imposed; e.g., conlab = 'All bar MS shocks inspl'
%-------------
% No output argument for this graph file
% See fn_seriesgraph.m, fn_forerrgraph.m.
%
% Tao Zha, March 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
vyrs = yfore(:,1);
hornum = cell(length(vyrs),1); % horizontal year (number)
count=0;
for k=vyrs'
count=count+1;
hornum{count}=num2str(k);
end
count=0;
for i = keyindx
count = count+1;
subplot(rnum,cnum,count)
plot(yacte(:,1)+yacte(:,2)/q_m,yacte(:,2+i),yfore(:,1)+yfore(:,2)/q_m,yfore(:,2+i),'--')
if (yfore(1,2)==0) % only for annual growth rates (not for, say, monthly annualized rates)
set(gca,'XLim',[vyrs(1) vyrs(end)])
set(gca,'XTick',vyrs)
set(gca,'XTickLabel',char(hornum))
end
if i==keyindx(1)
title(forelabel)
elseif i>=length(keyindx) %i>=length(keyindx)-1
xlabel(conlab)
end
%
grid
ylabel(char(ylab(i)))
end

60
matlab/ms-sbvar/cstz/fn_fprintmatrix.m
View File

@ -1,60 +0,0 @@
function fn_fprintmatrix(fid, M, nrows, ncols, indxFloat)
% Prints the matrix to an ascii file indexed by fid.
%
% Inputs:
% fid: Ascii file id. Example: fid = fopen('outdatainp_3s_stv_tvms6lags.prn','a');
% M: The matrix to be written to the file.
% nrows: Number of rows of M.
% ncols: Number of columns of M.
% indxFloat: 1 if double;
% 2 if single;
% 3 if only 3 significant digits
% 0 if integer.
%
% Copyright (C) 2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nrows~=size(M,1)
nrows
size(M,1)
error('fn_fprintmatrix(): Make sure the row number supplied match that of the matrix');
end
if ncols~=size(M,2)
ncols
size(M,2)
error('fn_fprintmatrix(): Make sure the column number supplied match that of the matrix');
end
for ki=1:nrows
for kj=1:ncols
if (indxFloat == 1)
fprintf(fid,' %.16e ',M((kj-1)*nrows+ki));
elseif (indxFloat == 2)
fprintf(fid,' %.8e ',M((kj-1)*nrows+ki));
elseif (indxFloat == 3)
fprintf(fid,' %.3e ',M((kj-1)*nrows+ki));
else
fprintf(fid,' %d ',M((kj-1)*nrows+ki));
end
if (kj==ncols)
fprintf(fid,'\n');
end
end
if (ki==nrows)
fprintf(fid,'\n\n');
end
end

58
matlab/ms-sbvar/cstz/fn_gfmean.m
View File

@ -1,58 +0,0 @@
function [Fmat,gvec] = fn_gfmean(b,P,Vi,nvar,ncoef,n0,np)
% [Fmat,gvec] = fn_gfmean(b,P,Vi,nvar,ncoef,n0,np)
%
% Mean of free lagged parameters g and original lagged parameters F, conditional on comtemporaneous b's
% See Waggoner and Zha's Gibbs sampling
%
% b: sum(n0)-element vector of mean estimate of A0 free parameters
% P: cell(nvar,1). In each cell, the transformation matrix that affects the posterior mean of A+ conditional on A0.
% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith
% equation lagged restriction matrix where k is a total of exogenous variables and
% ri is the number of free parameters. With this transformation, we have fi = Vi*gi
% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a
% vector of free parameters. There must be at least one free parameter left for
% the ith equation.
% nvar: number of endogeous variables
% ncoef: number of original lagged variables per equation
% n0: nvar-element vector, ith element represents the number of free A0 parameters in ith equation
% np: nvar-element vector, ith element represents the number of free A+ parameters in ith equation
%---------------
% Fmat: ncoef-by-nvar matrix of original lagged parameters A+. Column corresponding to equation.
% gvec: sum(np)-by-1 stacked vector of all free lagged parameters A+.
%
% Tao Zha, February 2000. Revised, August 2000.
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
b=b(:); n0=n0(:); np=np(:);
n0cum = [0;cumsum(n0)];
npcum = [0;cumsum(np)];
gvec = zeros(npcum(end),1);
Fmat = zeros(ncoef,nvar); % ncoef: maximum original lagged parameters per equation
if ~(length(b)==n0cum(end))
error('Make inputs n0 and length(b) match exactly')
end
for kj=1:nvar
bj = b(n0cum(kj)+1:n0cum(kj+1));
gj = P{kj}*bj;
gvec(npcum(kj)+1:npcum(kj+1)) = gj;
Fmat(:,kj) = Vi{kj}*gj;
end

83
matlab/ms-sbvar/cstz/fn_gibbsrvar.m
View File

@ -1,83 +0,0 @@
function [A0gbs, Wcell] = fn_gibbsrvar(A0gbs,UT,nvar,fss,n0,Indxcol)
% [A0gbs, Wcell] = fn_gibbsrvar(A0gbs,UT,nvar,fss,n0,Indxcol)
% One-step Gibbs sampler for restricted VARs in the structural form (including over-identified cases).
% Reference: "A Gibbs sampler for structural VARs" by D.F. Waggoner and T. Zha, ``
% Journal of Economic Dynamics & Control (JEDC) 28 (2003) 349-366.
% See Note Forecast (2) pp. 44-51, 70-71, and Theorem 1 and Section 3.1 in the WZ JEDC paper.
%
% A0gbs: the last draw of A0 matrix
% UT: cell(nvar,1) -- U_i*T_i in the proof of Theorem 1 where
% (1) a_i = U_i*b_i with b_i being a vector of free parameters
% (2) T_i (q_i-by-q_i) is from T_i*T_i'= S_i = inv(H0inv{i}/T) where H0inv is the inverse of
% the covariance martrix NOT divided by fss and S_i is defined in (14) on p.355 of the WZ JEDC paper.
% nvar: rank of A0 or # of variables
% fss: effective sample size == nSample (T)-lags+# of dummy observations
% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation
% Indxcol: a row vector indicating random columns this Gibbs draws.
% When this input is not supplied, the Gibbs draws all columns
%------------------
% A0bgs: nvar-by-nvar. New draw of A0 matrix in this Gibbs step
% Wcell: cell(nvar,1). In each cell, columns of Wcell{i} form an orthonormal basis w_1,...,w_qi in Section 3.1 in the WZ paper. Added 9/04.
%
% Written by Tao Zha, August 2000. Revised, September 2004.
% Copyright (c) by Waggoner and Zha
% Copyright (C) 2000-2011 Tao Zha and Daniel Waggoner
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if (nargin==5), Indxcol=[1:nvar]; end
%---------------- Local loop for Gibbs given last A0gbs ----------
Wcell = cell(length(Indxcol),1);
w = zeros(nvar,1);
for ki=Indxcol % given last A0gbs and generate new A0bgs
X = A0gbs; % WZ's Section 4.3
X(:,ki) = 0; % want to find non-zero sw s.t., X'*w=0
%*** Solving for w and getting an orthonormal basis. See Theorem 1 and also p.48 in Forecast II
[jL,Ux] = lu(X');
jIx0 = min(find(abs(diag(Ux))<eps)); % if isempty(jIx0), then something is wrong here
w(jIx0+1:end) = 0; % if jIx0+1>end, no effect on w0
w(jIx0) = 1;
jA = Ux(1:jIx0-1,1:jIx0-1);
jb = Ux(1:jIx0-1,jIx0);
jy = -jA\jb;
w(1:jIx0-1) = jy;
% Note: if jIx0=1 (which almost never happens for numerical stability reasons), no effect on w.
%*** Constructing orthonormal basis w_1, ..., w_qi at each Gibbs step
w0 = UT{ki}'*w;
w1 = w0/sqrt(sum(w0.^2));
[W,jnk] = qr(w1); % Columns of W form an orthonormal basis w_1,...,w_qi in Section 3.1 in the WZ paper
%*** Draw beta's according to Theorem 1 in the WZ JEDC paper.
gkbeta = zeros(n0(ki),1); % qi-by-1: greak beta's
jstd = sqrt(1/fss);
gkbeta(2:end) = jstd*randn(n0(ki)-1,1); % for beta_2, ..., beta_qi
%--- Unnormalized (i.e., not normalized) gamma or 1-d Wishart draw of beta_1 in Theorem 1 of the WZ JEDC paper.
jr = jstd*randn(fss+1,1);
if rand(1)<0.5
gkbeta(1) = sqrt(jr'*jr);
else
gkbeta(1) = -sqrt(jr'*jr);
end
%*** Getting a new ki_th column in A0
A0gbs(:,ki) = UT{ki}*(W*gkbeta); % n-by-1: U_i*T_i*W*beta;
Wcell{ki} = W; %q_i-by-1.
end

74
matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m
View File

@ -1,74 +0,0 @@
function [Tinv,UT,VHphalf,PU,VPU] = fn_gibbsrvar_setup(H0inv, Ui, Hpinv, Pmat, Vi, nvar, fss)
% [Tinv,UT,VHphalf,PU,VPU] = fn_gibbsrvar_setup.m(H0inv, Ui, Hpinv, Pmat, Vi, fss, nvar)
% Global setup outside the Gibbs loop to be used by fn_gibbsvar().
% Reference: "A Gibbs sampler for structural VARs" by D.F. Waggoner and T. Zha, ``
% Journal of Economic Dynamics & Control (JEDC) 28 (2003) 349-366.
% See Note Forecast (2) pp. 44-51, 70-71, and Theorem 1 and Section 3.1 in the WZ JEDC paper.
%
% H0inv: cell(nvar,1). Not divided by T yet. In each cell, inverse of posterior covariance matrix H0.
% The exponential term is b_i'*inv(H0)*b_i for the ith equation where b_i = U_i*a0_i.
% It resembles old SpH or Sbd in the exponent term in posterior of A0, but not divided by T yet.
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation. Imported from dnrprior.m.
% Hpinv: cell(nvar,1). In each cell, posterior inverse of covariance matrix Hp (A+) for the free parameters
% g_i = V_i*A+(:,i) in the ith equation.
% Pmat: cell(nvar,1). In each cell, the transformation matrix that affects the posterior mean of A+ conditional on A0.
% In other words, the posterior mean (of g_i) = Pmat{i}*b_i where g_i is a column vector of free parameters
% of A+(:,i)) given b_i (b_i is a column vector of free parameters of A0(:,i)).
% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith
% equation lagged restriction matrix where k (ncoef) is a total number of RHS variables and
% ri is the number of free parameters. With this transformation, we have fi = Vi*gi
% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a
% vector of free parameters. There must be at least one free parameter left for
% the ith equation. Imported from dnrprior.m.
% nvar: number of endogenous variables or rank of A0.
% fss: effective sample size (in the exponential term) = nSample - lags + ndobs (ndobs = # of dummy observations
% is set to 0 when fn_rnrprior_covres_dobs() is used where dummy observations are included as part of the explicit prior.
%-------------
% Tinv: cell(nvar,1). In each cell, inv(T_i) for T_iT_i'=S_i where S_i is defined on p.355 of the WZ JEDC paper.
% UT: cell(nvar,1). In each cell, U_i*T_i.
% VHphalf: cell(nvar,1). In each cell, V_i*sqrt(Hp_i).
% PU: cell(nvar,1). In each cell, Pmat{i}*U_i where Pmat{i} = P_i defined in (13) on p.353 of the WZ JEDC paper.
% VPU: cell(nvar,1). In each cell, V_i*P_i*U_i
%
% Written by Tao Zha, September 2004.
% Copyright (c) 2004 by Waggoner and Zha
% Copyright (C) 2004-2011 Tao Zha and Daniel Waggoner
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
%--- For A0.
Tinv = cell(nvar,1); % in each cell, inv(T_i) for T_iT_i'=S_i where S_i is defined on p.355 of the WZ JEDC paper.
UT = cell(nvar,1); % in each cell, U_i*T_i.
%--- For A+.
VHphalf = cell(nvar,1); % in each cell, V_i*sqrt(Hp_i).
PU = cell(nvar,1); % in each cell, Pmat{i}*U_i where Pmat{i} = P_i defined in (13) on p.353 of the WZ JEDC paper.
VPU = cell(nvar,1); % in each cell, V_i*P_i*U_i
%
for ki=1:nvar
%--- For A0.
Tinv{ki} = chol(H0inv{ki}/fss); % Tinv_i'*Tinv_i = inv(S_i) ==> T_i*T_i' = S_i where S_i = H0inv{i}/fss is defined on p.355 of the WZ JEDC paper.
UT{ki} = Ui{ki}/Tinv{ki}; % n-by-qi: U_i*T_i in (14) on p. 255 of the WZ JEDC paper.
%--- For A+.
VHphalf{ki} = Vi{ki}/chol(Hpinv{ki}); % where chol(Hpinv_i)*chol(Hpinv_i)'=Hpinv_i.
PU{ki} = Pmat{ki}*Ui{ki}';
VPU{ki} = Vi{ki}*PU{ki};
end

124
matlab/ms-sbvar/cstz/fn_imcgraph.m
View File

@ -1,124 +0,0 @@
function scaleout = fn_imcgraph(imf,nvar,imstp,xlab,ylab,indxGimfml,xTick)
% scaleout = fn_imcgraph(imf,nvar,imstp,xlab,ylab,indxGimfml,xTick)
% imcgraph: impulse, c (column: shock 1 to N), graph the ML point impulse response
%
% imf: imstp-by-nvar^2 matrix of impulse responses, column (responses to 1st shock, responses to 2nd shock
% etc), row (impusle steps),
% nvar: number of variables
% imstp: number of steps of impulse responses
% xlab,ylab: labels
% indxGimfml: 1, graph; 0, no graph
% xTick: optional. Eg: [12 24 36].
%---------------
% scaleout: column 1 represents maximums; column 2 minimums. Rows: nvar variables.
%
% NOTE: I added "indxGimfml" so this function may not be compatible with programs
% older than 03/06/99, TZ
%
% See imrgraph, fn_imcerrgraph, fn_imc2errgraph, imrerrgraph, fn_gyrfore in RVARcode
% Copyright (C) 1999-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin < 7, xTick = []; end
t = 1:imstp;
temp1=zeros(nvar,1);
temp2=zeros(nvar,1);
maxval=zeros(nvar,1);
minval=zeros(nvar,1);
for i = 1:nvar
for j = 1:nvar
temp1(j)=max(imf(:,(j-1)*nvar+i));
temp2(j)=min(imf(:,(j-1)*nvar+i));
end
maxval(i)=max(temp1);
minval(i)=min(temp2);
end
scaleout = [maxval(:) minval(:)];
%--------------
% Column j: Shock 1 to N; Row i: Responses to
%-------------
if indxGimfml
%figure
rowlabel = 1;
for i = 1:nvar % Responses of
columnlabel = 1;
if minval(i)<0
if maxval(i)<=0
yt=[minval(i) 0];
else
yt=[minval(i) 0 maxval(i)];
end
else % (minval(i) >=0)
if maxval(i) > 0
yt=[0 maxval(i)];
else % (identically zero responses)
yt=[-1 0 1];
end
end
scale=[1 imstp minval(i) maxval(i)];
for j = 1:nvar % To shocks
k1=(i-1)*nvar+j;
k2=(j-1)*nvar+i;
subplot(nvar,nvar,k1)
plot(t,imf(:,k2),t,zeros(length(imf(:,k2)),1),'r:');
set(gca,'XTick',xTick)
set(gca,'YTick',yt)
grid
axis(scale); % put limits on both axes.
%
% if maxval(i)>minval(i)
% set(gca,'YLim',[minval(i) maxval(i)])
% end
if isempty(xTick) %1 % No numbers on both axes
set(gca,'XTickLabel',' ');
set(gca,'YTickLabel',' ');
else % Put numbers on both axes
if i<nvar
set(gca,'XTickLabelMode','manual','XTickLabel',[])
end
if j>1
set(gca,'YTickLabel',' ');
end
end
if rowlabel == 1
%title(['x' num2str(j)])
%title(eval(['x' num2str(j)]))
title(char(xlab(j)))
end
if columnlabel == 1
%ylabel(['x' num2str(i)])
%ylabel(eval(['x' num2str(i)]))
ylabel(char(ylab(i)))
end
columnlabel = 0;
end
rowlabel = 0;
end
end

76
matlab/ms-sbvar/cstz/fn_impulse.m
View File

@ -1,76 +0,0 @@
function imf = fn_impulse(Bh,swish,nn)
% Computing impulse functions with
% imf = fn_impulse(Bh,swish,nn)
% imf is in a format that is the SAME as in RATS.
% Column: nvar responses to 1st shock,
% nvar responses to 2nd shock, and so on.
% Row: steps of impulse responses.
%-----------------
% Bh is the estimated reduced form coefficient in the form
% Y(T*nvar) = XB + U, X: T*k (may include all exogenous terms), B: k*nvar.
% The matrix form and dimension are the same as "Bh" from the function "sye.m";
% Column: 1st lag (with nvar variables) to lags (with nvar variables) + const = k.
% Note: columns correspond to equations.
% swish is the inv(A0) in the structural model y(t)A0 = e(t).
% Note: columns corresponding to equations.
% nn is the numbers of inputs [nvar,lags,# of steps of impulse responses].
%
% Written by Tao Zha.
% Copyright (c) 1994 by Tao Zha
% Copyright (C) 1994-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nvar = nn(1);
lags = nn(2);
imstep = nn(3); % number of steps for impulse responses
Ah = Bh';
% Row: nvar equations
% Column: 1st lag (with nvar variables) to lags (with nvar variables) + const = k.
imf = zeros(imstep,nvar*nvar);
% Column: nvar responses to 1st shock, nvar responses to 2nd shock, and so on.
% Row: steps of impulse responses.
M = zeros(nvar*(lags+1),nvar);
% Stack lags M's in the order of, e.g., [Mlags, ..., M2,M1;M0]
M(1:nvar,:) = swish';
Mtem = M(1:nvar,:); % temporary M.
% first (initial) responses to 1 standard deviation shock. Row: responses; Column: shocks
% * put in the form of "imf"
imf(1,:) = Mtem(:)';
t = 1;
ims1 = min([imstep-1 lags]);
while t <= ims1
Mtem = Ah(:,1:nvar*t)*M(1:nvar*t,:);
% Row: nvar equations, each for the nvar variables at tth lag
M(nvar+1:nvar*(t+1),:)=M(1:nvar*t,:);
M(1:nvar,:) = Mtem;
imf(t+1,:) = Mtem(:)';
% stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc.
t= t+1;
end
for t = lags+1:imstep-1
Mtem = Ah(:,1:nvar*lags)*M(1:nvar*lags,:);
% Row: nvar equations, each for the nvar variables at tth lag
M(nvar+1:nvar*(t+1),:) = M(1:nvar*t,:);
M(1:nvar,:)=Mtem;
imf(t+1,:) = Mtem(:)';
% stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc.
end

67
matlab/ms-sbvar/cstz/fn_rlrpostr.m
View File

@ -1,67 +0,0 @@
function [P,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0invtld,Hpinvtld,Ui,Vi)
% [P,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0tld,Hptld,Ui,Vi)
%
% Exporting random (i.e., random prior) Bayesian posterior matrices with linear restrictions
% See Waggoner and Zha's Gibbs sampling paper
%
% xtx: X'X: k-by-k where k=ncoef
% xty: X'Y: k-by-nvar
% yty: Y'Y: nvar-by-nvar
% Ptld: cell(nvar,1), transformation matrix that affects the (random walk) prior mean of A+ conditional on A0.
% H0invtld: cell(nvar,1), transformed inv covaraince for free parameters in A0(:,i).
% Hpinvtld: cell(nvar,1), transformed inv covaraince for free parameters in A+(:,i);
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation. Imported from dnrprior.m.
% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith
% equation lagged restriction matrix where k (ncoef) is a total number of RHS variables and
% ri is the number of free parameters. With this transformation, we have fi = Vi*gi
% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a
% vector of free parameters. There must be at least one free parameter left for
% the ith equation. Imported from dnrprior.m.
%-----------------
% P: cell(nvar,1). In each cell, the transformation matrix that affects the posterior mean of A+ conditional on A0.
% In other words, the posterior mean (of g_i) = P{i}*b_i where g_i is a column vector of free parameters
% of A+(:,i)) given b_i (b_i is a column vector of free parameters of A0(:,i)).
% H0inv: cell(nvar,1). Not divided by T yet. In each cell, inverse of posterior covariance matrix H0.
% The exponential term is b_i'*inv(H0)*b_i for the ith equation where b_i = U_i*a0_i.
% It resembles old SpH or Sbd in the exponent term in posterior of A0, but not divided by T yet.
% Hpinv: cell(nvar,1). In each cell, posterior inverse of covariance matrix Hp (A+) for the free parameters
% g_i = V_i*A+(:,i) in the ith equation.
%
% Tao Zha, February 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nvar = size(yty,1);
P = cell(nvar,1); % tld: tilda
H0inv = cell(nvar,1); % posterior inv(H0), resemble old SpH, but not divided by T yet.
Hpinv = cell(nvar,1); % posterior inv(Hp).
for n=1:nvar % one for each equation
Hpinv{n} = Vi{n}'*xtx*Vi{n} + Hpinvtld{n};
P1 = Vi{n}'*xty*Ui{n} + Hpinvtld{n}*Ptld{n};
P{n} = Hpinv{n}\P1;
H0inv{n} = Ui{n}'*yty*Ui{n} + H0invtld{n} + Ptld{n}'*Hpinvtld{n}*Ptld{n} ...
- P1'*P{n}; %P{n} = (Hpinv{n}\P1);
end

56
matlab/ms-sbvar/cstz/fn_rlrprior.m
View File

@ -1,56 +0,0 @@
function [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Ui,Vi,Pi,H0multi,Hpmulti,nvar)
% [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Ui,Vi,Pi,H0multi,Hpmulti,nvar)
%
% Exporting random Bayesian prior with linear restrictions
% See Waggoner and Zha's Gibbs sampling paper
%
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation. Imported from dnrprior.m.
% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith
% equation lagged restriction matrix where k (ncoef) is a total number of RHS variables and
% ri is the number of free parameters. With this transformation, we have fi = Vi*gi
% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a
% vector of free parameters. There must be at least one free parameter left for
% the ith equation. Imported from dnrprior.m.
% Pi: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations
% H0multi: nvar-by-nvar-by-nvar; H0 for different equations under asymmetric prior
% Hpmulti: ncoef-by-ncoef-by-nvar; H+ for different equations under asymmetric prior
% nvar: number of endogenous variables
% --------------------
% Ptld: cell(nvar,1). The prior mean of g_i is Ptld{i}*b_i;
% H0invtld: cell(nvar,1). Transformed inv covaraince for b_i, the free parameters in A0(:,i);
% Hpinvtld: cell(nvar,1). Transformed inv covaraince for g_i, the free parameters in A+(:,i);
%
% Tao Zha, February 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
Ptld = cell(nvar,1); % tld: tilda
H0invtld = cell(nvar,1); % H0 for different equations under linear restrictions
Hpinvtld = cell(nvar,1); % H+ for different equations under linear restrictions
for n=1:nvar % one for each equation
Hpinvtld{n} = Vi{n}'*(Hpmulti(:,:,n)\Vi{n});
Ptld{n} = (Hpinvtld{n}\Vi{n}')*(Hpmulti(:,:,n)\Pi)*Ui{n};
H0invtld{n} = Ui{n}'*(H0multi(:,:,n)\Ui{n}) + Ui{n}'*Pi'*(Hpmulti(:,:,n)\Pi)*Ui{n} ...
- Ptld{n}'*Hpinvtld{n}*Ptld{n};
end

297
matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m
View File

@ -1,297 +0,0 @@
function [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] ...
= fn_rnrprior_covres_dobs(nvar,q_m,lags,xdgel,mu,indxDummy,hpmsmd,indxmsmdeqn,nexo,asym0,asymp)
% Differs from fn_rnrprior_covres_dobs_tv(): no linear restrictions (Ui and Vi) have applied yet to this function, but
% linear restrictions are incorported in fn_rnrprior_covres_dobs_tv().
%
% Only works for the nexo=1 (constant term) case. To extend this to other exogenous variables, see fn_dataxy.m. 01/14/03.
% Differs from fn_rnrprior_covres.m in that dummy observations are included as part of the explicit prior. See Forcast II, pp.68-69b.
% More general than fn_rnrprior.m because when hpmsmd=0, fn_rnrprior_covres() is the same as fn_rnrprior().
% Allows for prior covariances for the MS and MD equations to achieve liquidity effects.
% Exports random Bayesian prior of Sims and Zha with asymmetric rior (but no linear restrictions yet)
% See Waggoner and Zha's Gibbs sampling paper and TVBVAR NOTES p. 71k.0.
%
% nvar: number of endogenous variables
% q_m: quarter or month
% lags: the maximum length of lag
% xdgel: T*nvar endogenous-variable matrix of raw or original data (no manipulation involved) with sample size including lags.
% Order of columns: (1) nvar endogenous variables; (2) constants will be automatically put in the last column.
% Used only to get variances of residuals for mu(1)-mu(5) and for dummy observations mu(5) and mu(6).
% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where
% mu(5) and mu(6) are NOT used here. See fn_dataxy.m for using mu(5) and mu(6).
% mu(1): overall tightness and also for A0; (0.57)
% mu(2): relative tightness for A+; (0.13)
% mu(3): relative tightness for the constant term; (0.1). NOTE: for other
% exogenous terms, the variance of each exogenous term must be taken into
% acount to eliminate the scaling factor.
% mu(4): tightness on lag decay; (1)
% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5)
% mu(6): weight on single dummy initial observation including constant
% (cointegration, unit roots, and stationarity); (5)
% indxDummy: 1: uses dummy observations to form part of an explicit prior; 0: no dummy observations as part of the prior.
% 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.
% 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.
% 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.
% This will give us a liquidity effect.
% indxmsmdeqn: 4-by-1 index for the locations of the MS and MD equation and for the locations of M and R.
% indxmsmdeqn(1) for MS and indxmsmdeqn(2) for MD.
% indxmsmdeqn(3) for M and indxmsmdeqn(4) for R.
% nexo: number of exogenous variables (if not specified, nexo=1 (constant) by default).
% The constant term is always put to the last of all endogenous and exogenous variables.
% asym0: nvar-by-nvar asymmetric prior on A0. Column -- equation.
% If ones(nvar,nvar), symmetric prior; if not, relative (asymmetric) tightness on A0.
% asymp: ncoef-1-by-nvar asymmetric prior on A+ bar constant. Column -- equation.
% If ones(ncoef-1,nvar), symmetric prior; if not, relative (asymmetric) tightness on A+.
% --------------------
% Pi: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations
% H0multi: nvar-by-nvar-by-nvar; H0 for different equations under asymmetric prior
% Hpmulti: ncoef-by-ncoef-by-nvar; H+ for different equations under asymmetric prior
% H0invmulti: nvar-by-nvar-by-nvar; inv(H0) for different equations under asymmetric prior
% Hpinvmulti: ncoef-by-ncoef-by-nvar; inv(H+) for different equations under asymmetric prior
%
% Tao Zha, February 2000. Revised, September 2000, February, May 2003.
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if (nargin<=8), nexo=1; end
ncoef = nvar*lags+nexo; % number of coefficients in *each* equation, RHS coefficients only.
H0multi=zeros(nvar,nvar,nvar); % H0 for different equations under asymmetric prior
Hpmulti=zeros(ncoef,ncoef,nvar); % H+ for different equations under asymmetric prior
H0invmulti=zeros(nvar,nvar,nvar); % inv(H0) for different equations under asymmetric prior
Hpinvmulti=zeros(ncoef,ncoef,nvar); % inv(H+) for different equations under asymmetric prior
%*** Constructing Pi for the ith equation under the random walk assumption
Pi = zeros(ncoef,nvar); % same for all equations
Pi(1:nvar,1:nvar) = eye(nvar); % random walk
%
%@@@ Prepared for Bayesian prior
%
%
% ** monthly lag decay in order to match quarterly decay: a*exp(bl) where
% ** l is the monthly lag. Suppose quarterly decay is 1/x where x=1,2,3,4.
% ** Let the decay of l1 (a*exp(b*l1)) match that of x1 (say, beginning: 1/1)
% ** and the decay of l2 (a*exp(b*l2)) match that of x2 (say, end: 1/5),
% ** we can solve for a and b which are
% ** b = (log_x1-log_x2)/(l1-l2), and a = x1*exp(-b*l1).
if q_m==12
l1 = 1; % 1st month == 1st quarter
xx1 = 1; % 1st quarter
l2 = lags; % last month
xx2 = 1/((ceil(lags/3))^mu(4)); % last quarter
%xx2 = 1/6; % last quarter
% 3rd quarter: i.e., we intend to let decay of the 6th month match
% that of the 3rd quarter, so that the 6th month decays a little
% faster than the second quarter which is 1/2.
if lags==1
b = 0;
else
b = (log(xx1)-log(xx2))/(l1-l2);
end
a = xx1*exp(-b*l1);
end
%
% *** specify the prior for each equation separately, SZ method,
% ** get the residuals from univariate regressions.
%
sgh = zeros(nvar,1); % square root
sgsh = sgh; % square
nSample=size(xdgel,1); % sample size-lags
yu = xdgel;
C = ones(nSample,1);
for k=1:nvar
[Bk,ek,junk1,junk2,junk3,junk4] = sye([yu(:,k) C],lags);
clear Bk junk1 junk2 junk3 junk4;
sgsh(k) = ek'*ek/(nSample-lags);
sgh(k) = sqrt(sgsh(k));
end
% ** prior variance for A0(:,1), same for all equations!!!
sg0bid = zeros(nvar,1); % Sigma0_bar diagonal only for the ith equation
for j=1:nvar
sg0bid(j) = 1/sgsh(j); % sgsh = sigmai^2
end
% ** prior variance for lagged and exogeous variables, same for all equations
sgpbid = zeros(ncoef,1); % Sigma_plus_bar, diagonal, for the ith equation
for i = 1:lags
if (q_m==12)
lagdecay = a*exp(b*i*mu(4));
end
%
for j = 1:nvar
if (q_m==12)
% exponential decay to match quarterly decay
sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation
elseif (q_m==4)
sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation
elseif (q_m==1)
sgpbid((i-1)*nvar+j) = (1/(i*4)^mu(4))^2/sgsh(j); % ith equation
else
error('Incompatibility with lags, check the possible errors!!!')
%warning('Incompatibility with lags, check the possible errors!!!')
%return
end
end
end
%
if indxDummy % Dummy observations as part of the explicit prior.
ndobs=nvar+1; % Number of dummy observations: nvar unit roots and 1 cointegration prior.
phibar = zeros(ndobs,ncoef);
%* constant term
const = ones(nvar+1,1);
const(1:nvar) = 0.0;
phibar(:,ncoef) = const; % the first nvar periods: no or zero constant!
xdgelint = mean(xdgel(1:lags,:),1); % mean of the first lags initial conditions
%* Dummies
for k=1:nvar
for m=1:lags
phibar(ndobs,nvar*(m-1)+k) = xdgelint(k);
phibar(k,nvar*(m-1)+k) = xdgelint(k);
% <<>> multiply hyperparameter later
end
end
phibar(1:nvar,:) = 1*mu(5)*phibar(1:nvar,:); % standard Sims and Zha prior
phibar(ndobs,:) = mu(6)*phibar(ndobs,:);
[phiq,phir]=qr(phibar,0);
xtxbar=phir'*phir; % phibar'*phibar. ncoef-by-ncoef. Reduced (not full) rank. See Forcast II, pp.69-69b.
end
%=================================================
% Computing the (prior) covariance matrix for A0, no data yet.
% As proved in pp.69a-69b, Forecast II, the following prior covariance of A0
% will remain the same after the dummy observations prior is incorporated.
% The dummy observation prior only affects the prior covariance of A+|A0.
% See pp.69a-69b for the proof.
%=================================================
%
%
% ** set up the conditional prior variance sg0bi and sgpbi.
sg0bida = mu(1)^2*sg0bid; % ith equation
sgpbida = mu(1)^2*mu(2)^2*sgpbid;
sgpbida(ncoef-nexo+1:ncoef) = mu(1)^2*mu(3)^2;
%<<>> No scaling adjustment has been made for exogenous terms other than constant
sgppbd = sgpbida(nvar+1:ncoef); % corresponding to A++, in a Sims-Zha paper
Hptd = zeros(ncoef);
Hptdi=Hptd;
Hptd(ncoef,ncoef)=sgppbd(ncoef-nvar);
Hptdinv(ncoef,ncoef)=1./sgppbd(ncoef-nvar);
% condtional on A0i, H_plus_tilde
if nargin<10 % the default is no asymmetric information
asym0 = ones(nvar,nvar); % if not ones, then we have relative (asymmetric) tightness
asymp = ones(ncoef-1,nvar); % for A+. Column -- equation
end
%**** Asymmetric Information
%asym0 = ones(nvar,nvar); % if not ones, then we have relative (asymmetric) tightness
%asymp = ones(ncoef-1,nvar); % pp: plus without constant. Column -- equation
%>>>>>> B: asymmetric prior variance for asymp <<<<<<<<
%
%for i = 1:lags
% rowif = (i-1)*nvar+1;
% rowil = i*nvar;
% idmatw0 = 0.5; % weight assigned to idmat0 in the formation of asymp
% if (i==1)
% asymp(rowif:rowil,:)=(1-idmatw0)*ones(nvar)+idmatw0*idmat0; % first lag
% % note: idmat1 is already transposed. Column -- equation
% else
% %asymp(rowif:rowil,1:nvar) = (1-idmatw0)*ones(nvar)+idmatw0*idmat0;
% % <<<<<<< toggle +
% % Note: already transposed, since idmat0 is transposed.
% % Meaning: column implies equation
% asymp(rowif:rowil,1:nvar) = ones(nvar);
% % >>>>>>> toggle -
% end
%end
%
%>>>>>> E: asymmetric prior variance for asymp <<<<<<<<
%=================================================
% Computing the final covariance matrix (S1,...,Sm) for the prior of A0,
% and final Gb=(G1,...,Gm) for A+ if asymmetric prior or for
% B if symmetric prior for A+
%=================================================
%
for i = 1:nvar
%------------------------------
% Introduce prior information on which variables "belong" in various equations.
% In this first trial, we just introduce this information here, in a model-specific way.
% Eventually this info has to be passed parametricly. In our first shot, we just damp down
% all coefficients except those on the diagonal.
%*** For A0
factor0=asym0(:,i);
sg0bd = sg0bida.*factor0; % Note, this only works for the prior variance Sg(i)
% of a0(i) being diagonal. If the prior variance Sg(i) is not
% diagonal, we have to the inverse to get inv(Sg(i)).
%sg0bdinv = 1./sg0bd;
% * unconditional variance on A0+
H0td = diag(sg0bd); % unconditional
%=== Correlation in the MS equation to get a liquidity effect.
if (i==indxmsmdeqn(1))
H0td(indxmsmdeqn(3),indxmsmdeqn(4)) = hpmsmd(1)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4)));
H0td(indxmsmdeqn(4),indxmsmdeqn(3)) = hpmsmd(1)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4)));
elseif (i==indxmsmdeqn(2))
H0td(indxmsmdeqn(3),indxmsmdeqn(4)) = hpmsmd(2)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4)));
H0td(indxmsmdeqn(4),indxmsmdeqn(3)) = hpmsmd(2)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4)));
end
H0tdinv = inv(H0td);
%H0tdinv = diag(sg0bdinv);
%
H0multi(:,:,i)=H0td;
H0invmulti(:,:,i)=H0tdinv;
%*** For A+
if ~(lags==0) % For A1 to remain random walk properties
factor1=asymp(1:nvar,i);
sg1bd = sgpbida(1:nvar).*factor1;
sg1bdinv = 1./sg1bd;
%
Hptd(1:nvar,1:nvar)=diag(sg1bd);
Hptdinv(1:nvar,1:nvar)=diag(sg1bdinv);
if lags>1
factorpp=asymp(nvar+1:ncoef-1,i);
sgpp_cbd = sgppbd(1:ncoef-nvar-1) .* factorpp;
sgpp_cbdinv = 1./sgpp_cbd;
Hptd(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbd);
Hptdinv(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbdinv);
% condtional on A0i, H_plus_tilde
end
end
%---------------
% The dummy observation prior affects only the prior covariance of A+|A0,
% but not the covariance of A0. See pp.69a-69b for the proof.
%---------------
if indxDummy % Dummy observations as part of the explicit prior.
Hpinvmulti(:,:,i)=Hptdinv + xtxbar;
Hpmulti(:,:,i) = inv(Hpinvmulti(:,:,i));
else
Hpmulti(:,:,i)=Hptd;
Hpinvmulti(:,:,i)=Hptdinv;
end
end

327
matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs_tv2.m
View File

@ -1,327 +0,0 @@
function [Pi_bar,H0tldcell_inv,Hptldcell_inv] ...
= fn_rnrprior_covres_dobs_tv2(nvar,nStates,indxScaleStates,q_m,lags,xdgel,mu,indxDummy,Ui,Vi,hpmsmd,indxmsmdeqn,nexo,asym0,asymp)
% Differs from fn_rnrprior_covres_dobs(): linear restrictions (Ui and Vi) have been incorported in fn_rnrprior_covres_dobs_tv?().
% Differs from fn_rnrprior_covres_dobs_tv(): allows an option to scale up the prior variance by nStates or not scale at all,
% so that the prior value is the same as the constant VAR when the parameters in all states are the same.
%
% Only works for the nexo=1 (constant term) case. To extend this to other exogenous variables, see fn_dataxy.m. 01/14/03.
% 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.
% Exports random Bayesian prior of Sims and Zha with asymmetric rior with linear restrictions already applied
% and with dummy observations (i.e., mu(5) and mu(6)) used as part of an explicit prior.
% This function allows for prior covariances for the MS and MD equations to achieve liquidity effects.
% See Waggoner and Zha's Gibbs sampling paper and TVBVAR NOTES pp. 71k.0 and 50-61.
%
% nvar: number of endogenous variables
% nStates: Number of states.
% indxScaleStates: if 0, no scale adjustment in the prior variance for the number of states in the function fn_rnrprior_covres_dobs_tv2();
% if 1: allows a scale adjustment, marking the prior variance bigger by the number of states.
% q_m: quarter or month
% lags: the maximum length of lag
% xdgel: T*nvar endogenous-variable matrix of raw or original data (no manipulation involved) with sample size including lags.
% Order of columns: (1) nvar endogenous variables; (2) constants will be automatically put in the last column.
% Used only to get variances of residuals for mu(1)-mu(5) and for dummy observations mu(5) and mu(6).
% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where
% mu(5) and mu(6) are NOT used here. See fn_dataxy.m for using mu(5) and mu(6).
% mu(1): overall tightness and also for A0; (0.57)
% mu(2): relative tightness for A+; (0.13)
% mu(3): relative tightness for the constant term; (0.1). NOTE: for other
% exogenous terms, the variance of each exogenous term must be taken into
% acount to eliminate the scaling factor.
% mu(4): tightness on lag decay; (1)
% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5)
% mu(6): weight on single dummy initial observation including constant
% (cointegration, unit roots, and stationarity); (5)
% NOTE: for this function, mu(5) and mu(6) are not used. See fn_dataxy.m for using mu(5) and mu(6).
% indxDummy: 1: uses dummy observations to form part of an explicit prior; 0: no dummy observations as part of the prior.
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi*si orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters
% within the state and si is the number of free states.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation in the order of [a_i for 1st state, ..., a_i for last state].
% Vi: nvar-by-1 cell. In each cell, k-by-ri*ti orthonormal basis for the null of the ith
% equation lagged restriction matrix where k is a total of exogenous variables and
% ri is the number of free parameters within the state and ti is the number of free states.
% With this transformation, we have fi = Vi*gi
% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a
% vector of free parameters. The ith equation is in the order of [nvar variables
% for 1st lag and 1st state, ..., nvar variables for last lag and 1st state, const for 1st state, nvar
% variables for 1st lag and 2nd state, nvar variables for last lag and 2nd state, const for 2nd state, and so on].
% 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.
% 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.
% 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.
% This will give us a liquidity effect. If hpmsmd=0, no such restrictions will be imposed.
% indxmsmdeqn: 4-by-1 index for the locations of the MS and MD equation and for the locations of M and R.
% indxmsmdeqn(1) for MS and indxmsmdeqn(2) for MD.
% indxmsmdeqn(3) for M and indxmsmdeqn(4) for R.
% nexo: number of exogenous variables (if not specified, nexo=1 (constant) by default).
% The constant term is always put to the last of all endogenous and exogenous variables.
% asym0: nvar-by-nvar asymmetric prior on A0. Column -- equation.
% If ones(nvar,nvar), symmetric prior; if not, relative (asymmetric) tightness on A0.
% asymp: ncoef-1-by-nvar asymmetric prior on A+ bar constant. Column -- equation.
% If ones(ncoef-1,nvar), symmetric prior; if not, relative (asymmetric) tightness on A+.
% --------------------
% Pi_bar: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations
% H0tldcell_inv: cell(nvar,1). The ith cell represents the ith equation, where the dim is
% qi*si-by-qi*si. The inverse of H0tld on p.60.
% Hptldcell_inv: cell(nvar,1). The ith cell represents the ith equation, where the dim is
% ri*ti-by-ri*ti.The inverse of Hptld on p.60.
%
% Differs from fn_rnrprior_covres_dobs(): linear restrictions (Ui and Vi) have been incorported in fn_rnrprior_covres_dobs_tv?().
% Differs from fn_rnrprior_covres_dobs_tv(): allows an option to scale up the prior variance by nStates or not scale at all.
% so that the prior value is the same as the constant VAR when the parameters in all states are the same.
% Tao Zha, February 2000. Revised, September 2000, 2001, February, May 2003, May 2004.
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<=12, nexo=1; end % <<>>1
ncoef = nvar*lags+nexo; % Number of coefficients in *each* equation for each state, RHS coefficients only.
ncoefsts = nStates*ncoef; % Number of coefficients in *each* equation in all states, RHS coefficients only.
H0tldcell_inv=cell(nvar,1); % inv(H0tilde) for different equations under asymmetric prior.
Hptldcell_inv=cell(nvar,1); % inv(H+tilde) for different equations under asymmetric prior.
%*** Constructing Pi_bar for the ith equation under the random walk assumption
Pi_bar = zeros(ncoef,nvar); % same for all equations
Pi_bar(1:nvar,1:nvar) = eye(nvar); % random walk
%
%@@@ Prepared for Bayesian prior
%
%
% ** monthly lag decay in order to match quarterly decay: a*exp(bl) where
% ** l is the monthly lag. Suppose quarterly decay is 1/x where x=1,2,3,4.
% ** Let the decay of l1 (a*exp(b*l1)) match that of x1 (say, beginning: 1/1)
% ** and the decay of l2 (a*exp(b*l2)) match that of x2 (say, end: 1/5),
% ** we can solve for a and b which are
% ** b = (log_x1-log_x2)/(l1-l2), and a = x1*exp(-b*l1).
if q_m==12
l1 = 1; % 1st month == 1st quarter
xx1 = 1; % 1st quarter
l2 = lags; % last month
xx2 = 1/((ceil(lags/3))^mu(4)); % last quarter
%xx2 = 1/6; % last quarter
% 3rd quarter: i.e., we intend to let decay of the 6th month match
% that of the 3rd quarter, so that the 6th month decays a little
% faster than the second quarter which is 1/2.
if lags==1
b = 0;
else
b = (log(xx1)-log(xx2))/(l1-l2);
end
a = xx1*exp(-b*l1);
end
%
% *** specify the prior for each equation separately, SZ method,
% ** get the residuals from univariate regressions.
%
sgh = zeros(nvar,1); % square root
sgsh = sgh; % square
nSample=size(xdgel,1); % sample size-lags
yu = xdgel;
C = ones(nSample,1);
for k=1:nvar
[Bk,ek,junk1,junk2,junk3,junk4] = sye([yu(:,k) C],lags);
clear Bk junk1 junk2 junk3 junk4;
sgsh(k) = ek'*ek/(nSample-lags);
sgh(k) = sqrt(sgsh(k));
end
% ** prior variance for A0(:,1), same for all equations!!!
sg0bid = zeros(nvar,1); % Sigma0_bar diagonal only for the ith equation
for j=1:nvar
sg0bid(j) = 1/sgsh(j); % sgsh = sigmai^2
end
% ** prior variance for lagged and exogeous variables, same for all equations
sgpbid = zeros(ncoef,1); % Sigma_plus_bar, diagonal, for the ith equation
for i = 1:lags
if (q_m==12)
lagdecay = a*exp(b*i*mu(4));
end
%
for j = 1:nvar
if (q_m==12)
% exponential decay to match quarterly decay
sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation
elseif (q_m==4)
sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation
elseif (q_m==1)
sgpbid((i-1)*nvar+j) = (1/(i*4)^mu(4))^2/sgsh(j); % ith equation
else
error('Incompatibility with lags, check the possible errors!!!')
%warning('Incompatibility with lags, check the possible errors!!!')
%return
end
end
end
%
if indxDummy % Dummy observations as part of the explicit prior.
ndobs=nvar+1; % Number of dummy observations: nvar unit roots and 1 cointegration prior.
phibar = zeros(ndobs,ncoef);
%* constant term
const = ones(nvar+1,1);
const(1:nvar) = 0.0;
phibar(:,ncoef) = const; % the first nvar periods: no or zero constant!
xdgelint = mean(xdgel(1:lags,:),1); % mean of the first lags initial conditions
%* Dummies
for k=1:nvar
for m=1:lags
phibar(ndobs,nvar*(m-1)+k) = xdgelint(k);
phibar(k,nvar*(m-1)+k) = xdgelint(k);
% <<>> multiply hyperparameter later
end
end
phibar(1:nvar,:) = 1*mu(5)*phibar(1:nvar,:); % standard Sims and Zha prior
phibar(ndobs,:) = mu(6)*phibar(ndobs,:);
[phiq,phir]=qr(phibar,0);
xtxbar=phir'*phir; % phibar'*phibar. ncoef-by-ncoef. Reduced (not full) rank. See Forcast II, pp.69-69b.
end
%=================================================
% Computing the (prior) covariance matrix for the posterior of A0, no data yet
%=================================================
%
%
% ** set up the conditional prior variance sg0bi and sgpbi.
sg0bida = mu(1)^2*sg0bid; % ith equation
sgpbida = mu(1)^2*mu(2)^2*sgpbid;
sgpbida(ncoef-nexo+1:ncoef) = mu(1)^2*mu(3)^2;
%<<>> No scaling adjustment has been made for exogenous terms other than constant
sgppbd = sgpbida(nvar+1:ncoef); % corresponding to A++, in a Sims-Zha paper
Hptd = zeros(ncoef);
Hptdi=Hptd;
Hptd(ncoef,ncoef)=sgppbd(ncoef-nvar);
Hptdinv(ncoef,ncoef)=1./sgppbd(ncoef-nvar);
% condtional on A0i, H_plus_tilde
if nargin<14 % <<>>1 Default is no asymmetric information
asym0 = ones(nvar,nvar); % if not ones, then we have relative (asymmetric) tightness
asymp = ones(ncoef-1,nvar); % for A+. Column -- equation
end
%**** Asymmetric Information
%asym0 = ones(nvar,nvar); % if not ones, then we have relative (asymmetric) tightness
%asymp = ones(ncoef-1,nvar); % pp: plus without constant. Column -- equation
%>>>>>> B: asymmetric prior variance for asymp <<<<<<<<
%
%for i = 1:lags
% rowif = (i-1)*nvar+1;
% rowil = i*nvar;
% idmatw0 = 0.5; % weight assigned to idmat0 in the formation of asymp
% if (i==1)
% asymp(rowif:rowil,:)=(1-idmatw0)*ones(nvar)+idmatw0*idmat0; % first lag
% % note: idmat1 is already transposed. Column -- equation
% else
% %asymp(rowif:rowil,1:nvar) = (1-idmatw0)*ones(nvar)+idmatw0*idmat0;
% % <<<<<<< toggle +
% % Note: already transposed, since idmat0 is transposed.
% % Meaning: column implies equation
% asymp(rowif:rowil,1:nvar) = ones(nvar);
% % >>>>>>> toggle -
% end
%end
%
%>>>>>> E: asymmetric prior variance for asymp <<<<<<<<
%=================================================
% Computing the final covariance matrix (S1,...,Sm) for the prior of A0,
% and final Gb=(G1,...,Gm) for A+ if asymmetric prior or for
% B if symmetric prior for A+
%=================================================
%
for i = 1:nvar
%------------------------------
% Introduce prior information on which variables "belong" in various equations.
% In this first trial, we just introduce this information here, in a model-specific way.
% Eventually this info has to be passed parametricly. In our first shot, we just damp down
% all coefficients except those on the diagonal.
%*** For A0
factor0=asym0(:,i);
sg0bd = sg0bida.*factor0; % Note, this only works for the prior variance Sg(i)
% of a0(i) being diagonal. If the prior variance Sg(i) is not
% diagonal, we have to the inverse to get inv(Sg(i)).
%sg0bdinv = 1./sg0bd;
% * unconditional variance on A0+
H0td = diag(sg0bd); % unconditional
%=== Correlation in the MS equation to get a liquidity effect.
if (i==indxmsmdeqn(1))
H0td(indxmsmdeqn(3),indxmsmdeqn(4)) = hpmsmd(1)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4)));
H0td(indxmsmdeqn(4),indxmsmdeqn(3)) = hpmsmd(1)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4)));
elseif (i==indxmsmdeqn(2))
H0td(indxmsmdeqn(3),indxmsmdeqn(4)) = hpmsmd(2)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4)));
H0td(indxmsmdeqn(4),indxmsmdeqn(3)) = hpmsmd(2)*sqrt(sg0bida(indxmsmdeqn(3))*sg0bida(indxmsmdeqn(4)));
end
H0tdinv = inv(H0td);
%H0tdinv = diag(sg0bdinv);
%
if indxScaleStates
H0tldcell_inv{i}=NaN;
else
H0tldcell_inv{i}=NaN;
end
%*** For A+
if ~(lags==0) % For A1 to remain random walk properties
factor1=asymp(1:nvar,i);
sg1bd = sgpbida(1:nvar).*factor1;
sg1bdinv = 1./sg1bd;
%
Hptd(1:nvar,1:nvar)=diag(sg1bd);
Hptdinv(1:nvar,1:nvar)=diag(sg1bdinv);
if lags>1
factorpp=asymp(nvar+1:ncoef-1,i);
sgpp_cbd = sgppbd(1:ncoef-nvar-1) .* factorpp;
sgpp_cbdinv = 1./sgpp_cbd;
Hptd(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbd);
Hptdinv(nvar+1:ncoef-1,nvar+1:ncoef-1)=diag(sgpp_cbdinv);
% condtional on A0i, H_plus_tilde
end
end
%---------------
% The dummy observation prior affects only the prior covariance of A+|A0,
% but not the covariance of A0. See pp.69a-69b for the proof.
%---------------
if indxDummy % Dummy observations as part of the explicit prior.
Hptdinv2 = Hptdinv + xtxbar; % Rename Hptdinv to Hptdinv2 because we want to keep Hptdinv diagonal in the next loop of i.
else
Hptdinv2 = Hptdinv;
end
if (indxScaleStates)
Hptldcell_inv{i}=NaN;
else
Hptldcell_inv{i}=NaN;
end
%Hptdinv_3 = kron(eye(nStates),Hptdinv); % ?????
end

42
matlab/ms-sbvar/cstz/fn_tran_a2b.m
View File

@ -1,42 +0,0 @@
function b = fn_tran_a2b(A0,Ui,nvar,n0)
% b = fn_tran_a2b(A0,Ui,nvar,n0)
% Transform A0 to free parameters b's. Note: columns correspond to equations
% See Waggoner and Zha's ``A Gibbs sampler for structural VARs''
%
% A0: nvar-by-nvar, contempareous matrix (columns correspond to equations)
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation.
% nvar: number of endogeous variables
% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation
%----------------
% b: sum(n0)-by-1 vector of A0 free parameters
%
% Tao Zha, February 2000. Revised, August 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
n0=n0(:);
n0cum = [0; cumsum(n0)];
b=zeros(n0cum(end),1);
for kj = 1:nvar
b(n0cum(kj)+1:n0cum(kj+1))=Ui{kj}'*A0(:,kj);
end

41
matlab/ms-sbvar/cstz/fn_tran_b2a.m
View File

@ -1,41 +0,0 @@
function A0 = fn_tran_b2a(b,Ui,nvar,n0)
% A0 = fn_tran_b2a(b,Ui,nvar,n0)
% Transform free parameters b's to A0. Note: columns correspond to equations
%
% b: sum(n0)-by-1 vector of A0 free parameters
% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith
% equation contemporaneous restriction matrix where qi is the number of free parameters.
% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector
% of total original parameters and bi is a vector of free parameters. When no
% restrictions are imposed, we have Ui = I. There must be at least one free
% parameter left for the ith equation.
% nvar: number of endogeous variables
% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation
%----------------
% A0: nvar-by-nvar, contempareous matrix (columns correspond to equations)
%
% Tao Zha, February 2000. Revised, August 2000.
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
b=b(:); n0=n0(:);
A0 = zeros(nvar);
n0cum = [0; cumsum(n0)];
for kj = 1:nvar
A0(:,kj) = Ui{kj}*b(n0cum(kj)+1:n0cum(kj+1));
end

46
matlab/ms-sbvar/cstz/fn_varoots.m
View File

@ -1,46 +0,0 @@
function rootsinv = fn_varoots(Bhat,nvar,lags)
%
% Using eigenvalues to find the inverse of all roots associated with the VAR proceess:
% y_t' = C + y_{t-1}'*B_1 + ... + Y_{t-p}'*B_p + u_t'.
% where columns correspond to equations. See also Judge (1), pp.753-755 where rows correspond to equations.
% Bhat: ncoef-by-nvar where ncoef=nvar*lags+nexo and nvar is the number of endogenous variables.
% Columns corresponds to equations with
% ncoef=[nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term]
% ..., nvar coef in the last lag, and nexo coefficients.
% Note that entries in the rows of Bhat that > nvar*lags are irrelevant.
% nvar: number of endogenous variables.
% lags: number of lags.
%-------
% rootsinv: a vector of nvar*lags inverse roots. When > 1, explosive. When all < 1, stationary.
%
% Tao Zha, September 2000
% Copyright (C) 2000-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if size(Bhat,1)<nvar*lags
disp(' ')
warning('Make sure that Bhat has at least nvar*lags rows')
return
end
%--------- Strack the VAR(p) to the VAR(1) with z_t = Az_{t-1}.
%
A1 = diag(ones(nvar*(lags-1),1));
A2 = [A1 zeros(nvar*(lags-1),nvar)];
A = [Bhat(1:nvar*lags,:)'; A2];
rootsinv=eig(A);

100
matlab/ms-sbvar/cstz/sye.m
View File

@ -1,100 +0,0 @@
function [Bh,e,xtx,xty,phi,y,ncoe,xr] = sye(z,lags)
% Now [Bh,e,xtx,xty,phi,y,ncoe,xr] = sye(z,lags)
% Old: [Bh,e,xtx,xty,phi,y,ncoe,Sigu,xtxinv] = sye(z,lags)
%
% Estimate a system of equations in the form of y(T*nvar) = XB + u,
% X--phi: T*k, B: k*nvar; where T=sp-lags, k=ncoe,
%
% z: (T+lags)-by-(nvar+1) raw data matrix (nvar of variables + constant).
% lags: number of lags
%--------------------
% Bh: k-by-nvar estimated reduced-form parameter; column: nvar;
% row: k=ncoe=[nvar for 1st lag, ..., nvar for last lag, const]
% e: estimated residual e = y -xBh, T-by-nvar
% xtx: X'X: k-by-k
% xty: X'Y: k-by-nvar
% phi: X; T-by-k; column: [nvar for 1st lag, ..., nvar for last lag, const]
% y: Y: T-by-nvar
% ncoe: number of coeffcients per equation: nvar*lags + 1
% xr: the economy size (k-by-k) in qr(phi) so that xr=chol(X'*X)
% Sigu: e'*e: nvar-by-nvar. Note, not divided (undivided) by "nobs"
% xtxinv: inv(X'X): k-by-k
%
% See also syed.m (allowing for more predetermined terms) which has not
% been yet updated as "sye.m".
%
% Note, "lags" is something I changed recently, so it may not be compatible
% with old programs, 10/15/98 by TAZ.
%
% Revised, 5/2/99. Replaced outputs Sigu and xtxinv with xr so that previous
% programs may be incompatible.
% Copyright (C) 1998-2011 Tao Zha
%
% This file is part of Dynare.
%
% Dynare is 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% ** setup of orders and lengths **
[sp,nvar] = size(z); % sp: sample period T include lags
nvar = nvar-1; % -1: takes out the counting of constant
ess = sp-lags; % effective sample size
sb = lags+1; % sample beginning
sl = sp; % sample last period
ncoe = nvar*lags + 1; % with constant
% ** construct X for Y = X*B + U where phi = X **
x = z(:,1:nvar);
C = z(:,nvar+1);
phi = zeros(ess,ncoe);
phi(:,ncoe) = C(1:ess);
for k=1:lags, phi(:,nvar*(k-1)+1:nvar*k) = x(sb-k:sl-k,:); end
% row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, const]
% Thus, # of columns is nvar*lags+1 = ncoef.
% ** estimate: B, XTX, residuals **
y = x(sb:sl,:);
%
%**** The following, though stable, is too slow *****
% [u d v]=svd(phi,0); %trial
% %xtx = phi'*phi; % X'X, k*k (ncoe*ncoe)
% vd=v.*(ones(size(v,2),1)*diag(d)'); %trial
% dinv = 1./diag(d); % inv(diag(d))
% vdinv=v.*(ones(size(v,2),1)*dinv'); %trial
% xtx=vd*vd';
% xtxinv = vdinv*vdinv';
% %xty = phi'*y; % X'Y
% uy = u'*y; %trial
% xty = vd*uy; %trial
% %Bh = xtx\xty; %inv(X'X)*(X'Y), k*m (ncoe*nvar).
% Bh = xtxinv*xty;
% %e = y - phi*Bh; % from Y = XB + U, e: (T-lags)*nvar
% e = y - u*uy;
%**** The following, though stable, is too slow *****
%===== (Fast but perhaps less accurate) alternative to the above =========
[xq,xr]=qr(phi,0);
xtx=xr'*xr;
xty=phi'*y;
Bh = xr\(xr'\xty);
e=y-phi*Bh;
%===== (Fast but perhaps less accurate) alternative to the above =========
%* not numerically stable way of computing e'*e
%Sigu = y'*y-xty'*Bh;
%Sigu = y'*(eye(ess)-phi*xtxinv*phi')*y; % probablly better way, commented out
% by TZ, 2/28/98. See following
%Sigu = y'*(eye(ess)-u*u')*y; % I think this is the best, TZ, 2/28/98
% Note, u*u'=x*inv(x'x)*x.