diff --git a/.gitmodules b/.gitmodules
index 9627637182f4da82a4b0d95de16639f9823190ca..91bcb4a1bb0585211366e9e15715f81f9785678f 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -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
diff --git a/contrib/ms-sbvar/TZcode b/contrib/ms-sbvar/TZcode
new file mode 160000
index 0000000000000000000000000000000000000000..1673c6f7170ac61ad3d3cb832eb3f038c8e4b3ef
--- /dev/null
+++ b/contrib/ms-sbvar/TZcode
@@ -0,0 +1 @@
+Subproject commit 1673c6f7170ac61ad3d3cb832eb3f038c8e4b3ef
diff --git a/license.txt b/license.txt
index d7ba91e074db67df4c66e8b0e892eb9363a55255..548610bf66400d8a325f1c9e5e594f171075c3d2 100644
--- a/license.txt
+++ b/license.txt
@@ -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+
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 6d0c90728f36584b973dbd4bed58b277be342d87..9eddd6ec2142792f0d979ec847bd5e407fe27288 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -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/'])
diff --git a/matlab/ms-sbvar/cstz/bfgsi.m b/matlab/ms-sbvar/cstz/bfgsi.m
deleted file mode 100644
index 3c2890e61ba09b1a883e4b5882d17428ccf0c6ed..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/bfgsi.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/csminit.m b/matlab/ms-sbvar/cstz/csminit.m
deleted file mode 100644
index 90b68239a08cb604181c60b959e28bbf6fd8b2ba..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/csminit.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/csminwel.m b/matlab/ms-sbvar/cstz/csminwel.m
deleted file mode 100644
index e3b7466d0ba0ec6503a8d11dce1f72275f8ebfac..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/csminwel.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_a0freefun.m b/matlab/ms-sbvar/cstz/fn_a0freefun.m
deleted file mode 100644
index 1ea2966a862e9018aeb94d5c61acf11475715108..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_a0freefun.m
+++ /dev/null
@@ -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;
diff --git a/matlab/ms-sbvar/cstz/fn_a0freegrad.m b/matlab/ms-sbvar/cstz/fn_a0freegrad.m
deleted file mode 100644
index 59e84c41833bbb258f6fb66b37ac04ccc0ef168b..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_a0freegrad.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_calyrqm.m b/matlab/ms-sbvar/cstz/fn_calyrqm.m
deleted file mode 100644
index a19d8169d8382e093e52ffeab93bee91d5aff552..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_calyrqm.m
+++ /dev/null
@@ -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);
diff --git a/matlab/ms-sbvar/cstz/fn_dataext.m b/matlab/ms-sbvar/cstz/fn_dataext.m
deleted file mode 100644
index 7df59f93ef108f62566eed8c3fd4925d50c3b78f..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_dataext.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_datana.m b/matlab/ms-sbvar/cstz/fn_datana.m
deleted file mode 100644
index 2da6d59a0ab06373348e46623e537cd1047caadf..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_datana.m
+++ /dev/null
@@ -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];
diff --git a/matlab/ms-sbvar/cstz/fn_dataxy.m b/matlab/ms-sbvar/cstz/fn_dataxy.m
deleted file mode 100644
index 164755312b0462bc1e9d937b6612268d040f187d..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_dataxy.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_ergodp.m b/matlab/ms-sbvar/cstz/fn_ergodp.m
deleted file mode 100644
index 7c5c70101f1255838aa2f94227ec76f3423cc59d..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_ergodp.m
+++ /dev/null
@@ -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));
diff --git a/matlab/ms-sbvar/cstz/fn_fcstidcnd.m b/matlab/ms-sbvar/cstz/fn_fcstidcnd.m
deleted file mode 100644
index 0328793d93fceeb79f89c6a047d0ce96a33ab88e..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_fcstidcnd.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_forecast.m b/matlab/ms-sbvar/cstz/fn_forecast.m
deleted file mode 100644
index 644131d7574ac22662d8c060f98c40d81138a217..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_forecast.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_foregraph.m b/matlab/ms-sbvar/cstz/fn_foregraph.m
deleted file mode 100644
index 4f4fc7ea6d18be6a0d7d3589bf99eeb9e1c32a49..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_foregraph.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_fprintmatrix.m b/matlab/ms-sbvar/cstz/fn_fprintmatrix.m
deleted file mode 100644
index 536cac6cc6fd3d0055b6e4d83b0d656454b90fb7..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_fprintmatrix.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_gfmean.m b/matlab/ms-sbvar/cstz/fn_gfmean.m
deleted file mode 100644
index f9c3dd4db0acc8a9a222361b5d531e26152ea1d2..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_gfmean.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_gibbsrvar.m b/matlab/ms-sbvar/cstz/fn_gibbsrvar.m
deleted file mode 100644
index cc5cd075c3996783e012017dbfb0ab5ad68210a3..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_gibbsrvar.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m b/matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m
deleted file mode 100644
index 54810c4a14cb8831a12ad20aafc45160b734046d..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_imcgraph.m b/matlab/ms-sbvar/cstz/fn_imcgraph.m
deleted file mode 100644
index 4eeebaaace83fd696fd7577099ba32f456aaccb6..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_imcgraph.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_impulse.m b/matlab/ms-sbvar/cstz/fn_impulse.m
deleted file mode 100644
index fb91120e97319d8c7b99d677c83eb6e5b06a3b78..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_impulse.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_rlrpostr.m b/matlab/ms-sbvar/cstz/fn_rlrpostr.m
deleted file mode 100644
index f8e5b12a5575402b0ddc336e6b9eb4721b35d6d9..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_rlrpostr.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_rlrprior.m b/matlab/ms-sbvar/cstz/fn_rlrprior.m
deleted file mode 100644
index 38600bee9f9a29af7005fadb524c82736d6fa02c..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_rlrprior.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m b/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m
deleted file mode 100644
index 41a22cab13132fd375817ed06fe8d335b8b58834..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m
+++ /dev/null
@@ -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
-
-
diff --git a/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs_tv2.m b/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs_tv2.m
deleted file mode 100644
index bc4fec59764aeac7eb5b383859143a6c7cf2d6f7..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs_tv2.m
+++ /dev/null
@@ -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
-
-
diff --git a/matlab/ms-sbvar/cstz/fn_tran_a2b.m b/matlab/ms-sbvar/cstz/fn_tran_a2b.m
deleted file mode 100644
index c10e2eee513a5d7163c4b19566f7d0b92776425c..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_tran_a2b.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_tran_b2a.m b/matlab/ms-sbvar/cstz/fn_tran_b2a.m
deleted file mode 100644
index 501b09326ad93c4d3702d106616ec790bb5b5fe9..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_tran_b2a.m
+++ /dev/null
@@ -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
diff --git a/matlab/ms-sbvar/cstz/fn_varoots.m b/matlab/ms-sbvar/cstz/fn_varoots.m
deleted file mode 100644
index 9020b216844e133d2a38c50ad94f75cb57592763..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/fn_varoots.m
+++ /dev/null
@@ -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);
diff --git a/matlab/ms-sbvar/cstz/sye.m b/matlab/ms-sbvar/cstz/sye.m
deleted file mode 100644
index b4a24326e4a6808800aeb8c6d34345ae9dfee180..0000000000000000000000000000000000000000
--- a/matlab/ms-sbvar/cstz/sye.m
+++ /dev/null
@@ -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.