From 6180bc2212ac00432b8a5fc9794e8c8a6a2dc188 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 19 Apr 2012 12:03:27 +0200 Subject: [PATCH] Add TZcode submodule to repository --- .gitmodules | 3 + contrib/ms-sbvar/TZcode | 1 + license.txt | 80 ++++- matlab/dynare_config.m | 2 +- matlab/ms-sbvar/cstz/bfgsi.m | 46 --- matlab/ms-sbvar/cstz/csminit.m | 216 ------------ matlab/ms-sbvar/cstz/csminwel.m | 306 ---------------- matlab/ms-sbvar/cstz/fn_a0freefun.m | 56 --- matlab/ms-sbvar/cstz/fn_a0freegrad.m | 59 ---- matlab/ms-sbvar/cstz/fn_calyrqm.m | 75 ---- matlab/ms-sbvar/cstz/fn_dataext.m | 60 ---- matlab/ms-sbvar/cstz/fn_datana.m | 117 ------- matlab/ms-sbvar/cstz/fn_dataxy.m | 164 --------- matlab/ms-sbvar/cstz/fn_ergodp.m | 33 -- matlab/ms-sbvar/cstz/fn_fcstidcnd.m | 322 ----------------- matlab/ms-sbvar/cstz/fn_forecast.m | 74 ---- matlab/ms-sbvar/cstz/fn_foregraph.m | 65 ---- matlab/ms-sbvar/cstz/fn_fprintmatrix.m | 60 ---- matlab/ms-sbvar/cstz/fn_gfmean.m | 58 ---- matlab/ms-sbvar/cstz/fn_gibbsrvar.m | 83 ----- matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m | 74 ---- matlab/ms-sbvar/cstz/fn_imcgraph.m | 124 ------- matlab/ms-sbvar/cstz/fn_impulse.m | 76 ---- matlab/ms-sbvar/cstz/fn_rlrpostr.m | 67 ---- matlab/ms-sbvar/cstz/fn_rlrprior.m | 56 --- .../ms-sbvar/cstz/fn_rnrprior_covres_dobs.m | 297 ---------------- .../cstz/fn_rnrprior_covres_dobs_tv2.m | 327 ------------------ matlab/ms-sbvar/cstz/fn_tran_a2b.m | 42 --- matlab/ms-sbvar/cstz/fn_tran_b2a.m | 41 --- matlab/ms-sbvar/cstz/fn_varoots.m | 46 --- matlab/ms-sbvar/cstz/sye.m | 100 ------ 31 files changed, 80 insertions(+), 3050 deletions(-) create mode 160000 contrib/ms-sbvar/TZcode delete mode 100644 matlab/ms-sbvar/cstz/bfgsi.m delete mode 100644 matlab/ms-sbvar/cstz/csminit.m delete mode 100644 matlab/ms-sbvar/cstz/csminwel.m delete mode 100644 matlab/ms-sbvar/cstz/fn_a0freefun.m delete mode 100644 matlab/ms-sbvar/cstz/fn_a0freegrad.m delete mode 100644 matlab/ms-sbvar/cstz/fn_calyrqm.m delete mode 100644 matlab/ms-sbvar/cstz/fn_dataext.m delete mode 100644 matlab/ms-sbvar/cstz/fn_datana.m delete mode 100644 matlab/ms-sbvar/cstz/fn_dataxy.m delete mode 100644 matlab/ms-sbvar/cstz/fn_ergodp.m delete mode 100644 matlab/ms-sbvar/cstz/fn_fcstidcnd.m delete mode 100644 matlab/ms-sbvar/cstz/fn_forecast.m delete mode 100644 matlab/ms-sbvar/cstz/fn_foregraph.m delete mode 100644 matlab/ms-sbvar/cstz/fn_fprintmatrix.m delete mode 100644 matlab/ms-sbvar/cstz/fn_gfmean.m delete mode 100644 matlab/ms-sbvar/cstz/fn_gibbsrvar.m delete mode 100644 matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m delete mode 100644 matlab/ms-sbvar/cstz/fn_imcgraph.m delete mode 100644 matlab/ms-sbvar/cstz/fn_impulse.m delete mode 100644 matlab/ms-sbvar/cstz/fn_rlrpostr.m delete mode 100644 matlab/ms-sbvar/cstz/fn_rlrprior.m delete mode 100644 matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m delete mode 100644 matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs_tv2.m delete mode 100644 matlab/ms-sbvar/cstz/fn_tran_a2b.m delete mode 100644 matlab/ms-sbvar/cstz/fn_tran_b2a.m delete mode 100644 matlab/ms-sbvar/cstz/fn_varoots.m delete mode 100644 matlab/ms-sbvar/cstz/sye.m diff --git a/.gitmodules b/.gitmodules index 962763718..91bcb4a1b 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 000000000..1673c6f71 --- /dev/null +++ b/contrib/ms-sbvar/TZcode @@ -0,0 +1 @@ +Subproject commit 1673c6f7170ac61ad3d3cb832eb3f038c8e4b3ef diff --git a/license.txt b/license.txt index d7ba91e07..548610bf6 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 6d0c90728..9eddd6ec2 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 3c2890e61..000000000 --- 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 . - -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 90b68239a..000000000 --- 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 . - -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 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 a1 - 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 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)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)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 e3b7466d0..000000000 --- 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 . - -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 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 1ea2966a8..000000000 --- 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 . - -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 59e84c418..000000000 --- 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 . - -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 a19d8169d..000000000 --- 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 . - -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 7df59f93e..000000000 --- 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 . - -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 2da6d59a0..000000000 --- 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 . - -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 164755312..000000000 --- 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 . - -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 7c5c70101..000000000 --- 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 . - -[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 0328793d9..000000000 --- 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 . - -DLSIdShock = ~isempty(eq_ms); % if not empty, the MS shock is identified as in DLS - -impsteps=size(imf3s_h,1); -if (forep # 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 - % 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 - % 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 644131d75..000000000 --- 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 . - -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 4f4fc7ea6..000000000 --- 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 . - -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 536cac6cc..000000000 --- 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 . - -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 f9c3dd4db..000000000 --- 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 . - -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 cc5cd075c..000000000 --- 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 . - -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))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 54810c4a1..000000000 --- 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 . - -%--- 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 4eeebaaac..000000000 --- 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 . - -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 i1 - 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 fb91120e9..000000000 --- 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 . - -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 f8e5b12a5..000000000 --- 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 . - -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 38600bee9..000000000 --- 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 . - -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 41a22cab1..000000000 --- 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. - -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 bc4fec597..000000000 --- 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. - -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 c10e2eee5..000000000 --- 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 . - -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 501b09326..000000000 --- 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 . - -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 9020b2168..000000000 --- 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 . - -if size(Bhat,1). - -% ** 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.