diff --git a/matlab/ms-sbvar/cstz/bfgsi.m b/matlab/ms-sbvar/cstz/bfgsi.m index d55a636281ac0b2c1a103e0a25da7f2e0cf6adee..3c2890e61ba09b1a883e4b5882d17428ccf0c6ed 100644 --- a/matlab/ms-sbvar/cstz/bfgsi.m +++ b/matlab/ms-sbvar/cstz/bfgsi.m @@ -1,46 +1,46 @@ -function H = bfgsi(H0,dg,dx) -% H = bfgsi(H0,dg,dx) -% dg is previous change in gradient; dx is previous change in x; -% 6/8/93 version that updates inverse hessian instead of hessian -% itself. -% Copyright by Christopher Sims 1996. This material may be freely -% reproduced and modified. -dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha) - -% Copyright (C) 1996-2011 Tao Zha and Christopher Sims -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if size(dg,2)>1 - dg=dg'; -end -if size(dx,2)>1 - dx=dx'; -end -Hdg = H0*dg; -dgdx = dg'*dx; -if (abs(dgdx) >1e-12) - H = H0 + (1+(dg'*Hdg)/dgdx)*(dx*dx')/dgdx - (dx*Hdg'+Hdg*dx')/dgdx; -else - if dispIndx - disp('bfgs update failed.') - disp(['|dg| = ' num2str(sqrt(dg'*dg)) '|dx| = ' num2str(sqrt(dx'*dx))]); - disp(['dg''*dx = ' num2str(dgdx)]) - disp(['|H*dg| = ' num2str(Hdg'*Hdg)]) - end - H=H0; -end -save H.dat H +function H = bfgsi(H0,dg,dx) +% H = bfgsi(H0,dg,dx) +% dg is previous change in gradient; dx is previous change in x; +% 6/8/93 version that updates inverse hessian instead of hessian +% itself. +% Copyright by Christopher Sims 1996. This material may be freely +% reproduced and modified. +dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha) + +% Copyright (C) 1996-2011 Tao Zha and Christopher Sims +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if size(dg,2)>1 + dg=dg'; +end +if size(dx,2)>1 + dx=dx'; +end +Hdg = H0*dg; +dgdx = dg'*dx; +if (abs(dgdx) >1e-12) + H = H0 + (1+(dg'*Hdg)/dgdx)*(dx*dx')/dgdx - (dx*Hdg'+Hdg*dx')/dgdx; +else + if dispIndx + disp('bfgs update failed.') + disp(['|dg| = ' num2str(sqrt(dg'*dg)) '|dx| = ' num2str(sqrt(dx'*dx))]); + disp(['dg''*dx = ' num2str(dgdx)]) + disp(['|H*dg| = ' num2str(Hdg'*Hdg)]) + end + H=H0; +end +save H.dat H diff --git a/matlab/ms-sbvar/cstz/csminit.m b/matlab/ms-sbvar/cstz/csminit.m index d1b52cbc62b17e8a95eab511fe9feb40d005f658..90b68239a08cb604181c60b959e28bbf6fd8b2ba 100644 --- a/matlab/ms-sbvar/cstz/csminit.m +++ b/matlab/ms-sbvar/cstz/csminit.m @@ -1,216 +1,216 @@ -function [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,varargin) -% [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,... -% P1,P2,P3,P4,P5,P6,P7,P8) -% retcodes: 0, normal step. 5, largest step still improves too fast. -% 4,2 back and forth adjustment of stepsize didn't finish. 3, smallest -% stepsize still improves too slow. 6, no improvement found. 1, zero -% gradient. -%--------------------- -% Modified 7/22/96 to omit variable-length P list, for efficiency and compilation. -% Places where the number of P's need to be altered or the code could be returned to -% its old form are marked with ARGLIST comments. -% -% Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs -% update. -% -% Fixed 7/19/93 to flip eigenvalues of H to get better performance when -% it's not psd. -% -% Fixed 02/19/05 to correct for low angle problems. -% -%tailstr = ')'; -%for i=nargin-6:-1:1 -% tailstr=[ ',P' num2str(i) tailstr]; -%end - -% Copyright (C) 1993-2011 Tao Zha and Christopher Sims -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha) - - -%ANGLE = .03; % when output of this variable becomes negative, we have wrong analytical graident -ANGLE = .005; % works for identified VARs and OLS -%THETA = .03; -THETA = .3; %(0<THETA<.5) THETA near .5 makes long line searches, possibly fewer iterations. -FCHANGE = 1000; -MINLAMB = 1e-9; -% fixed 7/15/94 -% MINDX = .0001; -% MINDX = 1e-6; -MINDFAC = .01; -fcount=0; -lambda=1; -xhat=x0; -f=f0; -fhat=f0; -g = g0; -gnorm = norm(g); -% -if (gnorm < 1.e-12) && ~badg % put ~badg 8/4/94 - retcode =1; - dxnorm=0; - % gradient convergence -else - % with badg true, we don't try to match rate of improvement to directional - % derivative. We're satisfied just to get some improvement in f. - % - %if(badg) - % dx = -g*FCHANGE/(gnorm*gnorm); - % dxnorm = norm(dx); - % if dxnorm > 1e12 - % disp('Bad, small gradient problem.') - % dx = dx*FCHANGE/dxnorm; - % end - %else - % Gauss-Newton step; - %---------- Start of 7/19/93 mod --------------- - %[v d] = eig(H0); - %toc - %d=max(1e-10,abs(diag(d))); - %d=abs(diag(d)); - %dx = -(v.*(ones(size(v,1),1)*d'))*(v'*g); -% toc - dx = -H0*g; -% toc - dxnorm = norm(dx); - if dxnorm > 1e12 - if dispIndx, disp('Near-singular H problem.'), end - dx = dx*FCHANGE/dxnorm; - end - dfhat = dx'*g0; - %end - % - % - if ~badg - % test for alignment of dx with gradient and fix if necessary - a = -dfhat/(gnorm*dxnorm); - if a<ANGLE - dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g; - % suggested alternate code: --------------------- - dx = dx*dxnorm/norm(dx); % Added 02/19/05 by CAS. This keeps scale invariant to the angle correction - % ------------------------------------------------ - dfhat = dx'*g; - % dxnorm = norm(dx); % Removed 02/19/05 by CAS. This line unnecessary with modification that keeps scale invariant - if dispIndx, disp(sprintf('Correct for low angle: %g',a)), end - end - end - if dispIndx, disp(sprintf('Predicted improvement: %18.9f',-dfhat/2)), end - % - % Have OK dx, now adjust length of step (lambda) until min and - % max improvement rate criteria are met. - done=0; - factor=3; - shrink=1; - lambdaMin=0; - lambdaMax=inf; - lambdaPeak=0; - fPeak=f0; - lambdahat=0; - while ~done - if size(x0,2)>1 - dxtest=x0+dx'*lambda; - else - dxtest=x0+dx*lambda; - end - % home - f = feval(fcn,dxtest,varargin{:}); - %ARGLIST - %f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); - % f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8); - if dispIndx, disp(sprintf('lambda = %10.5g; f = %20.7f',lambda,f )), end - %debug - %disp(sprintf('Improvement too great? f0-f: %g, criterion: %g',f0-f,-(1-THETA)*dfhat*lambda)) - if f<fhat - fhat=f; - xhat=dxtest; - lambdahat = lambda; - end - fcount=fcount+1; - shrinkSignal = (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | (badg & (f0-f) < 0) ; - growSignal = ~badg & ( (lambda > 0) & (f0-f > -(1-THETA)*dfhat*lambda) ); - if shrinkSignal && ( (lambda>lambdaPeak) || (lambda<0) ) - if (lambda>0) && ((~shrink) || (lambda/factor <= lambdaPeak)) - shrink=1; - factor=factor^.6; - while lambda/factor <= lambdaPeak - factor=factor^.6; - end - %if (abs(lambda)*(factor-1)*dxnorm < MINDX) || (abs(lambda)*(factor-1) < MINLAMB) - if abs(factor-1)<MINDFAC - if abs(lambda)<4 - retcode=2; - else - retcode=7; - end - done=1; - end - end - if (lambda<lambdaMax) && (lambda>lambdaPeak) - lambdaMax=lambda; - end - lambda=lambda/factor; - if abs(lambda) < MINLAMB - if (lambda > 0) && (f0 <= fhat) - % try going against gradient, which may be inaccurate - if dispIndx, lambda = -lambda*factor^6, end - else - if lambda < 0 - retcode = 6; - else - retcode = 3; - end - done = 1; - end - end - elseif (growSignal && lambda>0) || (shrinkSignal && ((lambda <= lambdaPeak) && (lambda>0))) - if shrink - shrink=0; - factor = factor^.6; - %if ( abs(lambda)*(factor-1)*dxnorm< MINDX ) || ( abs(lambda)*(factor-1)< MINLAMB) - if abs(factor-1)<MINDFAC - if abs(lambda)<4 - retcode=4; - else - retcode=7; - end - done=1; - end - end - if ( f<fPeak ) && (lambda>0) - fPeak=f; - lambdaPeak=lambda; - if lambdaMax<=lambdaPeak - lambdaMax=lambdaPeak*factor*factor; - end - end - lambda=lambda*factor; - if abs(lambda) > 1e20; - retcode = 5; - done =1; - end - else - done=1; - if factor < 1.2 - retcode=7; - else - retcode=0; - end - end - end -end -if dispIndx, disp(sprintf('Norm of dx %10.5g', dxnorm)), end +function [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,varargin) +% [fhat,xhat,fcount,retcode] = csminit(fcn,x0,f0,g0,badg,H0,... +% P1,P2,P3,P4,P5,P6,P7,P8) +% retcodes: 0, normal step. 5, largest step still improves too fast. +% 4,2 back and forth adjustment of stepsize didn't finish. 3, smallest +% stepsize still improves too slow. 6, no improvement found. 1, zero +% gradient. +%--------------------- +% Modified 7/22/96 to omit variable-length P list, for efficiency and compilation. +% Places where the number of P's need to be altered or the code could be returned to +% its old form are marked with ARGLIST comments. +% +% Fixed 7/17/93 to use inverse-hessian instead of hessian itself in bfgs +% update. +% +% Fixed 7/19/93 to flip eigenvalues of H to get better performance when +% it's not psd. +% +% Fixed 02/19/05 to correct for low angle problems. +% +%tailstr = ')'; +%for i=nargin-6:-1:1 +% tailstr=[ ',P' num2str(i) tailstr]; +%end + +% Copyright (C) 1993-2011 Tao Zha and Christopher Sims +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha) + + +%ANGLE = .03; % when output of this variable becomes negative, we have wrong analytical graident +ANGLE = .005; % works for identified VARs and OLS +%THETA = .03; +THETA = .3; %(0<THETA<.5) THETA near .5 makes long line searches, possibly fewer iterations. +FCHANGE = 1000; +MINLAMB = 1e-9; +% fixed 7/15/94 +% MINDX = .0001; +% MINDX = 1e-6; +MINDFAC = .01; +fcount=0; +lambda=1; +xhat=x0; +f=f0; +fhat=f0; +g = g0; +gnorm = norm(g); +% +if (gnorm < 1.e-12) && ~badg % put ~badg 8/4/94 + retcode =1; + dxnorm=0; + % gradient convergence +else + % with badg true, we don't try to match rate of improvement to directional + % derivative. We're satisfied just to get some improvement in f. + % + %if(badg) + % dx = -g*FCHANGE/(gnorm*gnorm); + % dxnorm = norm(dx); + % if dxnorm > 1e12 + % disp('Bad, small gradient problem.') + % dx = dx*FCHANGE/dxnorm; + % end + %else + % Gauss-Newton step; + %---------- Start of 7/19/93 mod --------------- + %[v d] = eig(H0); + %toc + %d=max(1e-10,abs(diag(d))); + %d=abs(diag(d)); + %dx = -(v.*(ones(size(v,1),1)*d'))*(v'*g); +% toc + dx = -H0*g; +% toc + dxnorm = norm(dx); + if dxnorm > 1e12 + if dispIndx, disp('Near-singular H problem.'), end + dx = dx*FCHANGE/dxnorm; + end + dfhat = dx'*g0; + %end + % + % + if ~badg + % test for alignment of dx with gradient and fix if necessary + a = -dfhat/(gnorm*dxnorm); + if a<ANGLE + dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g; + % suggested alternate code: --------------------- + dx = dx*dxnorm/norm(dx); % Added 02/19/05 by CAS. This keeps scale invariant to the angle correction + % ------------------------------------------------ + dfhat = dx'*g; + % dxnorm = norm(dx); % Removed 02/19/05 by CAS. This line unnecessary with modification that keeps scale invariant + if dispIndx, disp(sprintf('Correct for low angle: %g',a)), end + end + end + if dispIndx, disp(sprintf('Predicted improvement: %18.9f',-dfhat/2)), end + % + % Have OK dx, now adjust length of step (lambda) until min and + % max improvement rate criteria are met. + done=0; + factor=3; + shrink=1; + lambdaMin=0; + lambdaMax=inf; + lambdaPeak=0; + fPeak=f0; + lambdahat=0; + while ~done + if size(x0,2)>1 + dxtest=x0+dx'*lambda; + else + dxtest=x0+dx*lambda; + end + % home + f = feval(fcn,dxtest,varargin{:}); + %ARGLIST + %f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); + % f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8); + if dispIndx, disp(sprintf('lambda = %10.5g; f = %20.7f',lambda,f )), end + %debug + %disp(sprintf('Improvement too great? f0-f: %g, criterion: %g',f0-f,-(1-THETA)*dfhat*lambda)) + if f<fhat + fhat=f; + xhat=dxtest; + lambdahat = lambda; + end + fcount=fcount+1; + shrinkSignal = (~badg & (f0-f < max([-THETA*dfhat*lambda 0]))) | (badg & (f0-f) < 0) ; + growSignal = ~badg & ( (lambda > 0) & (f0-f > -(1-THETA)*dfhat*lambda) ); + if shrinkSignal && ( (lambda>lambdaPeak) || (lambda<0) ) + if (lambda>0) && ((~shrink) || (lambda/factor <= lambdaPeak)) + shrink=1; + factor=factor^.6; + while lambda/factor <= lambdaPeak + factor=factor^.6; + end + %if (abs(lambda)*(factor-1)*dxnorm < MINDX) || (abs(lambda)*(factor-1) < MINLAMB) + if abs(factor-1)<MINDFAC + if abs(lambda)<4 + retcode=2; + else + retcode=7; + end + done=1; + end + end + if (lambda<lambdaMax) && (lambda>lambdaPeak) + lambdaMax=lambda; + end + lambda=lambda/factor; + if abs(lambda) < MINLAMB + if (lambda > 0) && (f0 <= fhat) + % try going against gradient, which may be inaccurate + if dispIndx, lambda = -lambda*factor^6, end + else + if lambda < 0 + retcode = 6; + else + retcode = 3; + end + done = 1; + end + end + elseif (growSignal && lambda>0) || (shrinkSignal && ((lambda <= lambdaPeak) && (lambda>0))) + if shrink + shrink=0; + factor = factor^.6; + %if ( abs(lambda)*(factor-1)*dxnorm< MINDX ) || ( abs(lambda)*(factor-1)< MINLAMB) + if abs(factor-1)<MINDFAC + if abs(lambda)<4 + retcode=4; + else + retcode=7; + end + done=1; + end + end + if ( f<fPeak ) && (lambda>0) + fPeak=f; + lambdaPeak=lambda; + if lambdaMax<=lambdaPeak + lambdaMax=lambdaPeak*factor*factor; + end + end + lambda=lambda*factor; + if abs(lambda) > 1e20; + retcode = 5; + done =1; + end + else + done=1; + if factor < 1.2 + retcode=7; + else + retcode=0; + end + end + end +end +if dispIndx, disp(sprintf('Norm of dx %10.5g', dxnorm)), end diff --git a/matlab/ms-sbvar/cstz/csminwel.m b/matlab/ms-sbvar/cstz/csminwel.m index f335ad8a56adf59a85ba9b8ee251d87a74275bd5..e3b7466d0ba0ec6503a8d11dce1f72275f8ebfac 100644 --- a/matlab/ms-sbvar/cstz/csminwel.m +++ b/matlab/ms-sbvar/cstz/csminwel.m @@ -1,306 +1,306 @@ -function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,varargin) -%[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit,varargin) -% fcn: string naming the objective function to be minimized -% x0: initial value of the parameter vector -% H0: initial value for the inverse Hessian. Must be positive definite. -% grad: Either a string naming a function that calculates the gradient, or the null matrix. -% If it's null, the program calculates a numerical gradient. In this case fcn must -% be written so that it can take a matrix argument and produce a row vector of values. -% crit: Convergence criterion. Iteration will cease when it proves impossible to improve the -% function value by more than crit. -% nit: Maximum number of iterations. -% varargin: A list of optional length of additional parameters that get handed off to fcn each -% time it is called. -% Note that if the program ends abnormally, it is possible to retrieve the current x, -% f, and H from the files g1.mat and H.mat that are written at each iteration and at each -% hessian update, respectively. (When the routine hits certain kinds of difficulty, it -% write g2.mat and g3.mat as well. If all were written at about the same time, any of them -% may be a decent starting point. One can also start from the one with best function value.) -% NOTE: The display on screen can be turned off by seeting dispIndx=0 in this -% function. This option is used for the loop operation. T. Zha, 2 May 2000 -% NOTE: You may want to change stps to 1.0e-02 or 1.0e-03 to get a better convergence. August, 2006 - -% Copyright (C) 1993-2011 Tao Zha and Christopher Sims -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -Verbose = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha) -dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha) - -[nx,no]=size(x0); -nx=max(nx,no); -NumGrad= ( ~ischar(grad) | length(grad)==0); -done=0; -itct=0; -fcount=0; -snit=100; -%tailstr = ')'; -%stailstr = []; -% Lines below make the number of Pi's optional. This is inefficient, though, and precludes -% use of the matlab compiler. Without them, we use feval and the number of Pi's must be -% changed with the editor for each application. Places where this is required are marked -% with ARGLIST comments -%for i=nargin-6:-1:1 -% tailstr=[ ',P' num2str(i) tailstr]; -% stailstr=[' P' num2str(i) stailstr]; -%end -f0 = eval([fcn '(x0,varargin{:})']); -%ARGLIST -%f0 = feval(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); -% disp('first fcn in csminwel.m ----------------') % Jinill on 9/5/95 -if f0 > 1e50, disp('Bad initial parameter.'), return, end -if NumGrad - if length(grad)==0 - [g badg] = numgradcd(fcn,x0, varargin{:}); - %ARGLIST - %[g badg] = numgradcd(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); - else - badg=any(find(grad==0)); - g=grad; - end - %numgradcd(fcn,x0,P1,P2,P3,P4); -else - [g badg] = eval([grad '(x0,varargin{:})']); - %ARGLIST - %[g badg] = feval(grad,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); -end -retcode3=101; -x=x0; -f=f0; -H=H0; -cliff=0; -while ~done - g1=[]; g2=[]; g3=[]; - %addition fj. 7/6/94 for control - if dispIndx - disp('-----------------') - disp('-----------------') - %disp('f and x at the beginning of new iteration') - disp(sprintf('f at the beginning of new iteration, %20.10f',f)) - %-----------Comment out this line if the x vector is long---------------- - disp([sprintf('x = ') sprintf('%15.8g%15.8g%15.8g%15.8g%15.8g\n',x)]); - end - %------------------------- - itct=itct+1; - [f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,varargin{:}); - %ARGLIST - %[f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,P1,P2,P3,P4,P5,P6,P7,... - % P8,P9,P10,P11,P12,P13); - % itct=itct+1; - fcount = fcount+fc; - % erased on 8/4/94 - % if (retcode == 1) || (abs(f1-f) < crit) - % done=1; - % end - % if itct > nit - % done = 1; - % retcode = -retcode; - % end - if retcode1 ~= 1 - if retcode1==2 || retcode1==4 - wall1=1; badg1=1; - else - if NumGrad - [g1 badg1] = numgradcd(fcn, x1,varargin{:}); - %ARGLIST - %[g1 badg1] = numgradcd(fcn, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,... - % P10,P11,P12,P13); - else - [g1 badg1] = eval([grad '(x1,varargin{:})']); - %ARGLIST - %[g1 badg1] = feval(grad, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,... - % P10,P11,P12,P13); - end - wall1=badg1; - % g1 - save g1.mat g1 x1 f1 varargin; - %ARGLIST - %save g1 g1 x1 f1 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13; - end - if wall1 % && (~done) by Jinill - % Bad gradient or back and forth on step length. Possibly at - % cliff edge. Try perturbing search direction. - % - %fcliff=fh;xcliff=xh; - if dispIndx - disp(' ') - disp('************************* Random search. *****************************************') - disp('************************* Random search. *****************************************') - disp(' ') - pause(1.0) - end - Hcliff=H+diag(diag(H).*rand(nx,1)); - if dispIndx, disp('Cliff. Perturbing search direction.'), end - [f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,varargin{:}); - %ARGLIST - %[f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,P1,P2,P3,P4,... - % P5,P6,P7,P8,P9,P10,P11,P12,P13); - fcount = fcount+fc; % put by Jinill - if f2 < f - if retcode2==2 || retcode2==4 - wall2=1; badg2=1; - else - if NumGrad - [g2 badg2] = numgradcd(fcn, x2,varargin{:}); - %ARGLIST - %[g2 badg2] = numgradcd(fcn, x2,P1,P2,P3,P4,P5,P6,P7,P8,... - % P9,P10,P11,P12,P13); - else - [g2 badg2] = eval([grad '(x2,varargin{:})']); - %ARGLIST - %[g2 badg2] = feval(grad,x2,P1,P2,P3,P4,P5,P6,P7,P8,... - % P9,P10,P11,P12,P13); - end - wall2=badg2; - % g2 - if dispIndx, badg2, end - save g2.mat g2 x2 f2 varargin - %ARGLIST - %save g2 g2 x2 f2 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13; - end - if wall2 - if dispIndx, disp('Cliff again. Try traversing'), end - if norm(x2-x1) < 1e-13 - f3=f; x3=x; badg3=1;retcode3=101; - else - gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1); - [f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),varargin{:}); - %ARGLIST - %[f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),P1,P2,P3,... - % P4,P5,P6,P7,P8,... - % P9,P10,P11,P12,P13); - fcount = fcount+fc; % put by Jinill - if retcode3==2 || retcode3==4 - wall3=1; badg3=1; - else - if NumGrad - [g3 badg3] = numgradcd(fcn, x3,varargin{:}); - %ARGLIST - %[g3 badg3] = numgradcd(fcn, x3,P1,P2,P3,P4,P5,P6,P7,P8,... - % P9,P10,P11,P12,P13); - else - [g3 badg3] = eval([grad '(x3,varargin{:})']); - %ARGLIST - %[g3 badg3] = feval(grad,x3,P1,P2,P3,P4,P5,P6,P7,P8,... - % P9,P10,P11,P12,P13); - end - wall3=badg3; - % g3 - if dispIndx, badg3, end - save g3.mat g3 x3 f3 varargin; - %ARGLIST - %save g3 g3 x3 f3 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13; - end - end - else - f3=f; x3=x; badg3=1; retcode3=101; - end - else - f3=f; x3=x; badg3=1;retcode3=101; - end - else - % normal iteration, no walls, or else we're finished here. - f2=f; f3=f; badg2=1; badg3=1; retcode2=101; retcode3=101; - end - else - f1=f; f2=f; f3=f; retcode2=retcode1; retcode3=retcode1; - end - %how to pick gh and xh - if f3<f && badg3==0 - if dispIndx, ih=3, end - fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3; - elseif f2<f && badg2==0 - if dispIndx, ih=2, end - fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2; - elseif f1<f && badg1==0 - if dispIndx, ih=1, end - fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1; - else - [fh,ih] = min([f1,f2,f3]); - if dispIndx, disp(sprintf('ih = %d',ih)), end - %eval(['xh=x' num2str(ih) ';']) - switch ih - case 1 - xh=x1; - case 2 - xh=x2; - case 3 - xh=x3; - end %case - %eval(['gh=g' num2str(ih) ';']) - %eval(['retcodeh=retcode' num2str(ih) ';']) - retcodei=[retcode1,retcode2,retcode3]; - retcodeh=retcodei(ih); - if exist('gh') - nogh=isempty(gh); - else - nogh=1; - end - if nogh - if NumGrad - [gh badgh] = feval('numgrad',fcn,xh,varargin{:}); - else - [gh badgh] = feval('grad', xh,varargin{:}); - end - end - badgh=1; - end - %end of picking - %ih - %fh - %xh - %gh - %badgh - stuck = (abs(fh-f) < crit); - if (~badg)&(~badgh)&(~stuck) - H = bfgsi(H,gh-g,xh-x); - end - if Verbose - if dispIndx - disp('----') - disp(sprintf('Improvement on iteration %d = %18.9f',itct,f-fh)) - end - end - - if itct > nit - if dispIndx, disp('iteration count termination'), end - done = 1; - elseif stuck - if dispIndx, disp('improvement < crit termination'), end - done = 1; - end - rc=retcodeh; - if rc == 1 - if dispIndx, disp('zero gradient'), end - elseif rc == 6 - if dispIndx, disp('smallest step still improving too slow, reversed gradient'), end - elseif rc == 5 - if dispIndx, disp('largest step still improving too fast'), end - elseif (rc == 4) || (rc==2) - if dispIndx, disp('back and forth on step length never finished'), end - elseif rc == 3 - if dispIndx, disp('smallest step still improving too slow'), end - elseif rc == 7 - if dispIndx, disp('warning: possible inaccuracy in H matrix'), end - end - - f=fh; - x=xh; - g=gh; - badg=badgh; -end -% what about making an m-file of 10 lines including numgrad.m -% since it appears three times in csminwel.m +function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,varargin) +%[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit,varargin) +% fcn: string naming the objective function to be minimized +% x0: initial value of the parameter vector +% H0: initial value for the inverse Hessian. Must be positive definite. +% grad: Either a string naming a function that calculates the gradient, or the null matrix. +% If it's null, the program calculates a numerical gradient. In this case fcn must +% be written so that it can take a matrix argument and produce a row vector of values. +% crit: Convergence criterion. Iteration will cease when it proves impossible to improve the +% function value by more than crit. +% nit: Maximum number of iterations. +% varargin: A list of optional length of additional parameters that get handed off to fcn each +% time it is called. +% Note that if the program ends abnormally, it is possible to retrieve the current x, +% f, and H from the files g1.mat and H.mat that are written at each iteration and at each +% hessian update, respectively. (When the routine hits certain kinds of difficulty, it +% write g2.mat and g3.mat as well. If all were written at about the same time, any of them +% may be a decent starting point. One can also start from the one with best function value.) +% NOTE: The display on screen can be turned off by seeting dispIndx=0 in this +% function. This option is used for the loop operation. T. Zha, 2 May 2000 +% NOTE: You may want to change stps to 1.0e-02 or 1.0e-03 to get a better convergence. August, 2006 + +% Copyright (C) 1993-2011 Tao Zha and Christopher Sims +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +Verbose = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha) +dispIndx = 0; % 1: turn on all the diplays on the screen; 0: turn off (Added by T. Zha) + +[nx,no]=size(x0); +nx=max(nx,no); +NumGrad= ( ~ischar(grad) | length(grad)==0); +done=0; +itct=0; +fcount=0; +snit=100; +%tailstr = ')'; +%stailstr = []; +% Lines below make the number of Pi's optional. This is inefficient, though, and precludes +% use of the matlab compiler. Without them, we use feval and the number of Pi's must be +% changed with the editor for each application. Places where this is required are marked +% with ARGLIST comments +%for i=nargin-6:-1:1 +% tailstr=[ ',P' num2str(i) tailstr]; +% stailstr=[' P' num2str(i) stailstr]; +%end +f0 = eval([fcn '(x0,varargin{:})']); +%ARGLIST +%f0 = feval(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); +% disp('first fcn in csminwel.m ----------------') % Jinill on 9/5/95 +if f0 > 1e50, disp('Bad initial parameter.'), return, end +if NumGrad + if length(grad)==0 + [g badg] = numgradcd(fcn,x0, varargin{:}); + %ARGLIST + %[g badg] = numgradcd(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); + else + badg=any(find(grad==0)); + g=grad; + end + %numgradcd(fcn,x0,P1,P2,P3,P4); +else + [g badg] = eval([grad '(x0,varargin{:})']); + %ARGLIST + %[g badg] = feval(grad,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); +end +retcode3=101; +x=x0; +f=f0; +H=H0; +cliff=0; +while ~done + g1=[]; g2=[]; g3=[]; + %addition fj. 7/6/94 for control + if dispIndx + disp('-----------------') + disp('-----------------') + %disp('f and x at the beginning of new iteration') + disp(sprintf('f at the beginning of new iteration, %20.10f',f)) + %-----------Comment out this line if the x vector is long---------------- + disp([sprintf('x = ') sprintf('%15.8g%15.8g%15.8g%15.8g%15.8g\n',x)]); + end + %------------------------- + itct=itct+1; + [f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,varargin{:}); + %ARGLIST + %[f1 x1 fc retcode1] = csminit(fcn,x,f,g,badg,H,P1,P2,P3,P4,P5,P6,P7,... + % P8,P9,P10,P11,P12,P13); + % itct=itct+1; + fcount = fcount+fc; + % erased on 8/4/94 + % if (retcode == 1) || (abs(f1-f) < crit) + % done=1; + % end + % if itct > nit + % done = 1; + % retcode = -retcode; + % end + if retcode1 ~= 1 + if retcode1==2 || retcode1==4 + wall1=1; badg1=1; + else + if NumGrad + [g1 badg1] = numgradcd(fcn, x1,varargin{:}); + %ARGLIST + %[g1 badg1] = numgradcd(fcn, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,... + % P10,P11,P12,P13); + else + [g1 badg1] = eval([grad '(x1,varargin{:})']); + %ARGLIST + %[g1 badg1] = feval(grad, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,... + % P10,P11,P12,P13); + end + wall1=badg1; + % g1 + save g1.mat g1 x1 f1 varargin; + %ARGLIST + %save g1 g1 x1 f1 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13; + end + if wall1 % && (~done) by Jinill + % Bad gradient or back and forth on step length. Possibly at + % cliff edge. Try perturbing search direction. + % + %fcliff=fh;xcliff=xh; + if dispIndx + disp(' ') + disp('************************* Random search. *****************************************') + disp('************************* Random search. *****************************************') + disp(' ') + pause(1.0) + end + Hcliff=H+diag(diag(H).*rand(nx,1)); + if dispIndx, disp('Cliff. Perturbing search direction.'), end + [f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,varargin{:}); + %ARGLIST + %[f2 x2 fc retcode2] = csminit(fcn,x,f,g,badg,Hcliff,P1,P2,P3,P4,... + % P5,P6,P7,P8,P9,P10,P11,P12,P13); + fcount = fcount+fc; % put by Jinill + if f2 < f + if retcode2==2 || retcode2==4 + wall2=1; badg2=1; + else + if NumGrad + [g2 badg2] = numgradcd(fcn, x2,varargin{:}); + %ARGLIST + %[g2 badg2] = numgradcd(fcn, x2,P1,P2,P3,P4,P5,P6,P7,P8,... + % P9,P10,P11,P12,P13); + else + [g2 badg2] = eval([grad '(x2,varargin{:})']); + %ARGLIST + %[g2 badg2] = feval(grad,x2,P1,P2,P3,P4,P5,P6,P7,P8,... + % P9,P10,P11,P12,P13); + end + wall2=badg2; + % g2 + if dispIndx, badg2, end + save g2.mat g2 x2 f2 varargin + %ARGLIST + %save g2 g2 x2 f2 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13; + end + if wall2 + if dispIndx, disp('Cliff again. Try traversing'), end + if norm(x2-x1) < 1e-13 + f3=f; x3=x; badg3=1;retcode3=101; + else + gcliff=((f2-f1)/((norm(x2-x1))^2))*(x2-x1); + [f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),varargin{:}); + %ARGLIST + %[f3 x3 fc retcode3] = csminit(fcn,x,f,gcliff,0,eye(nx),P1,P2,P3,... + % P4,P5,P6,P7,P8,... + % P9,P10,P11,P12,P13); + fcount = fcount+fc; % put by Jinill + if retcode3==2 || retcode3==4 + wall3=1; badg3=1; + else + if NumGrad + [g3 badg3] = numgradcd(fcn, x3,varargin{:}); + %ARGLIST + %[g3 badg3] = numgradcd(fcn, x3,P1,P2,P3,P4,P5,P6,P7,P8,... + % P9,P10,P11,P12,P13); + else + [g3 badg3] = eval([grad '(x3,varargin{:})']); + %ARGLIST + %[g3 badg3] = feval(grad,x3,P1,P2,P3,P4,P5,P6,P7,P8,... + % P9,P10,P11,P12,P13); + end + wall3=badg3; + % g3 + if dispIndx, badg3, end + save g3.mat g3 x3 f3 varargin; + %ARGLIST + %save g3 g3 x3 f3 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13; + end + end + else + f3=f; x3=x; badg3=1; retcode3=101; + end + else + f3=f; x3=x; badg3=1;retcode3=101; + end + else + % normal iteration, no walls, or else we're finished here. + f2=f; f3=f; badg2=1; badg3=1; retcode2=101; retcode3=101; + end + else + f1=f; f2=f; f3=f; retcode2=retcode1; retcode3=retcode1; + end + %how to pick gh and xh + if f3<f && badg3==0 + if dispIndx, ih=3, end + fh=f3;xh=x3;gh=g3;badgh=badg3;retcodeh=retcode3; + elseif f2<f && badg2==0 + if dispIndx, ih=2, end + fh=f2;xh=x2;gh=g2;badgh=badg2;retcodeh=retcode2; + elseif f1<f && badg1==0 + if dispIndx, ih=1, end + fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1; + else + [fh,ih] = min([f1,f2,f3]); + if dispIndx, disp(sprintf('ih = %d',ih)), end + %eval(['xh=x' num2str(ih) ';']) + switch ih + case 1 + xh=x1; + case 2 + xh=x2; + case 3 + xh=x3; + end %case + %eval(['gh=g' num2str(ih) ';']) + %eval(['retcodeh=retcode' num2str(ih) ';']) + retcodei=[retcode1,retcode2,retcode3]; + retcodeh=retcodei(ih); + if exist('gh') + nogh=isempty(gh); + else + nogh=1; + end + if nogh + if NumGrad + [gh badgh] = feval('numgrad',fcn,xh,varargin{:}); + else + [gh badgh] = feval('grad', xh,varargin{:}); + end + end + badgh=1; + end + %end of picking + %ih + %fh + %xh + %gh + %badgh + stuck = (abs(fh-f) < crit); + if (~badg)&(~badgh)&(~stuck) + H = bfgsi(H,gh-g,xh-x); + end + if Verbose + if dispIndx + disp('----') + disp(sprintf('Improvement on iteration %d = %18.9f',itct,f-fh)) + end + end + + if itct > nit + if dispIndx, disp('iteration count termination'), end + done = 1; + elseif stuck + if dispIndx, disp('improvement < crit termination'), end + done = 1; + end + rc=retcodeh; + if rc == 1 + if dispIndx, disp('zero gradient'), end + elseif rc == 6 + if dispIndx, disp('smallest step still improving too slow, reversed gradient'), end + elseif rc == 5 + if dispIndx, disp('largest step still improving too fast'), end + elseif (rc == 4) || (rc==2) + if dispIndx, disp('back and forth on step length never finished'), end + elseif rc == 3 + if dispIndx, disp('smallest step still improving too slow'), end + elseif rc == 7 + if dispIndx, disp('warning: possible inaccuracy in H matrix'), end + end + + f=fh; + x=xh; + g=gh; + badg=badgh; +end +% what about making an m-file of 10 lines including numgrad.m +% since it appears three times in csminwel.m diff --git a/matlab/ms-sbvar/cstz/fn_a0freefun.m b/matlab/ms-sbvar/cstz/fn_a0freefun.m index 0f800783478f83a46cc0ef301b92efaf7dec0e59..1ea2966a862e9018aeb94d5c61acf11475715108 100644 --- a/matlab/ms-sbvar/cstz/fn_a0freefun.m +++ b/matlab/ms-sbvar/cstz/fn_a0freefun.m @@ -1,56 +1,56 @@ -function of = fn_a0freefun(b,Ui,nvar,n0,fss,H0inv) -% of = fn_a0freefun(b,Ui,nvar,n0,fss,H0inv) -% -% Negative logPosterior function for squeesed A0 free parameters, which are b's in the WZ notation -% Note: columns correspond to equations -% -% b: sum(n0)-by-1 vector of A0 free parameters -% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith -% equation contemporaneous restriction matrix where qi is the number of free parameters. -% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector -% of total original parameters and bi is a vector of free parameters. When no -% restrictions are imposed, we have Ui = I. There must be at least one free -% parameter left for the ith equation. -% nvar: number of endogeous variables -% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation -% fss: nSample-lags (plus ndobs if dummies are included) -% H0inv: cell(nvar,1). In each cell, posterior inverse of covariance inv(H0) for the ith equation, -% resembling old SpH in the exponent term in posterior of A0, but not divided by T yet. -%---------------- -% of: objective function (negative logPosterior) -% -% Tao Zha, February 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -b=b(:); n0=n0(:); - -A0 = zeros(nvar); -n0cum = [0;cumsum(n0)]; -tra = 0.0; -for kj = 1:nvar - bj = b(n0cum(kj)+1:n0cum(kj+1)); - A0(:,kj) = Ui{kj}*bj; - tra = tra + 0.5*bj'*H0inv{kj}*bj; % negative exponential term -end - -[A0l,A0u] = lu(A0); - -ada = -fss*sum(log(abs(diag(A0u)))); % negative log determinant of A0 raised to power T - -of = ada + tra; +function of = fn_a0freefun(b,Ui,nvar,n0,fss,H0inv) +% of = fn_a0freefun(b,Ui,nvar,n0,fss,H0inv) +% +% Negative logPosterior function for squeesed A0 free parameters, which are b's in the WZ notation +% Note: columns correspond to equations +% +% b: sum(n0)-by-1 vector of A0 free parameters +% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith +% equation contemporaneous restriction matrix where qi is the number of free parameters. +% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector +% of total original parameters and bi is a vector of free parameters. When no +% restrictions are imposed, we have Ui = I. There must be at least one free +% parameter left for the ith equation. +% nvar: number of endogeous variables +% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation +% fss: nSample-lags (plus ndobs if dummies are included) +% H0inv: cell(nvar,1). In each cell, posterior inverse of covariance inv(H0) for the ith equation, +% resembling old SpH in the exponent term in posterior of A0, but not divided by T yet. +%---------------- +% of: objective function (negative logPosterior) +% +% Tao Zha, February 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +b=b(:); n0=n0(:); + +A0 = zeros(nvar); +n0cum = [0;cumsum(n0)]; +tra = 0.0; +for kj = 1:nvar + bj = b(n0cum(kj)+1:n0cum(kj+1)); + A0(:,kj) = Ui{kj}*bj; + tra = tra + 0.5*bj'*H0inv{kj}*bj; % negative exponential term +end + +[A0l,A0u] = lu(A0); + +ada = -fss*sum(log(abs(diag(A0u)))); % negative log determinant of A0 raised to power T + +of = ada + tra; diff --git a/matlab/ms-sbvar/cstz/fn_a0freegrad.m b/matlab/ms-sbvar/cstz/fn_a0freegrad.m index 9ad0e205215230c0d47786c15cc3028e78556227..59e84c41833bbb258f6fb66b37ac04ccc0ef168b 100644 --- a/matlab/ms-sbvar/cstz/fn_a0freegrad.m +++ b/matlab/ms-sbvar/cstz/fn_a0freegrad.m @@ -1,59 +1,59 @@ -function [g,badg] = fn_a0freegrad(b,Ui,nvar,n0,fss,H0inv) -% [g,badg] = a0freegrad(b,Ui,nvar,n0,fss,H0inv) -% Analytical gradient for a0freefun.m in use of csminwel.m. See Dhrymes's book. -% -% b: sum(n0)-by-1 vector of A0 free parameters -% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith -% equation contemporaneous restriction matrix where qi is the number of free parameters. -% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector -% of total original parameters and bi is a vector of free parameters. When no -% restrictions are imposed, we have Ui = I. There must be at least one free -% parameter left for the ith equation. -% nvar: number of endogeous variables -% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation -% fss: nSample-lags (plus ndobs if dummies are included) -% H0inv: cell(nvar,1). In each cell, posterior inverse of covariance inv(H0) for the ith equation, -% resembling old SpH in the exponent term in posterior of A0, but not divided by T yet. -%--------------- -% g: sum(n0)-by-1 analytical gradient for a0freefun.m -% badg: 0, the value that is used in csminwel.m -% -% Tao Zha, February 2000. Revised, August 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -b=b(:); n0 = n0(:); - -A0 = zeros(nvar); -n0cum = [0;cumsum(n0)]; -g = zeros(n0cum(end),1); -badg = 0; - -%*** The derivative of the exponential term w.r.t. each free paramater -for kj = 1:nvar - bj = b(n0cum(kj)+1:n0cum(kj+1)); - g(n0cum(kj)+1:n0cum(kj+1)) = H0inv{kj}*bj; - A0(:,kj) = Ui{kj}*bj; -end -B=inv(A0'); - -%*** Add the derivative of -Tlog|A0| w.r.t. each free paramater -for ki = 1:sum(n0) - n = max(find( (ki-n0cum)>0 )); % note, 1<=n<=nvar - g(ki) = g(ki) - fss*B(:,n)'*Ui{n}(:,ki-n0cum(n)); -end +function [g,badg] = fn_a0freegrad(b,Ui,nvar,n0,fss,H0inv) +% [g,badg] = a0freegrad(b,Ui,nvar,n0,fss,H0inv) +% Analytical gradient for a0freefun.m in use of csminwel.m. See Dhrymes's book. +% +% b: sum(n0)-by-1 vector of A0 free parameters +% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith +% equation contemporaneous restriction matrix where qi is the number of free parameters. +% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector +% of total original parameters and bi is a vector of free parameters. When no +% restrictions are imposed, we have Ui = I. There must be at least one free +% parameter left for the ith equation. +% nvar: number of endogeous variables +% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation +% fss: nSample-lags (plus ndobs if dummies are included) +% H0inv: cell(nvar,1). In each cell, posterior inverse of covariance inv(H0) for the ith equation, +% resembling old SpH in the exponent term in posterior of A0, but not divided by T yet. +%--------------- +% g: sum(n0)-by-1 analytical gradient for a0freefun.m +% badg: 0, the value that is used in csminwel.m +% +% Tao Zha, February 2000. Revised, August 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +b=b(:); n0 = n0(:); + +A0 = zeros(nvar); +n0cum = [0;cumsum(n0)]; +g = zeros(n0cum(end),1); +badg = 0; + +%*** The derivative of the exponential term w.r.t. each free paramater +for kj = 1:nvar + bj = b(n0cum(kj)+1:n0cum(kj+1)); + g(n0cum(kj)+1:n0cum(kj+1)) = H0inv{kj}*bj; + A0(:,kj) = Ui{kj}*bj; +end +B=inv(A0'); + +%*** Add the derivative of -Tlog|A0| w.r.t. each free paramater +for ki = 1:sum(n0) + n = max(find( (ki-n0cum)>0 )); % note, 1<=n<=nvar + g(ki) = g(ki) - fss*B(:,n)'*Ui{n}(:,ki-n0cum(n)); +end diff --git a/matlab/ms-sbvar/cstz/fn_calyrqm.m b/matlab/ms-sbvar/cstz/fn_calyrqm.m index e8d07fd8b32f93f51d72151fbfa792079b17bb6d..a19d8169d8382e093e52ffeab93bee91d5aff552 100644 --- a/matlab/ms-sbvar/cstz/fn_calyrqm.m +++ b/matlab/ms-sbvar/cstz/fn_calyrqm.m @@ -1,75 +1,75 @@ -function [Myrqm,nMyrqm] = fn_calyrqm(q_m,Byrqm,Eyrqm) -% [Myrqm,nMyrqm] = fn_calyrqm(q_m,Byrqm,Eyrqm) -% -% Given the beginning and end years and quarters (months), export a matrix of all years and -% quarters (months) for these years and in between -% -% q_m: 4 if quarterly and 12 if monthly -% Byrqm: [year quarter(month)] -- all integers, the begining year and quarter (month) -% Eyrqm: [year quarter(month)] -- all integers, the end year and quarter (month) -%------------------- -% Myrqm: matrix of all years and quarters (months) between and incl. Byrqm and Eyrqm -% nMyrqm: number of data points incl. Byrqm and Eyrqm -% -% Tao Zha, April 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if ~isempty(find(Byrqm-round(Byrqm))) | (q_m-round(q_m)) | ~isempty(find(Byrqm-round(Byrqm))) - error('argin qm, Byrqm, or Eyrqm must of integer') -elseif Byrqm(1)>Eyrqm(1) - error('Eyrqm(1) must be equal to or greater than Byrqm(1)') -elseif Byrqm(1)==Eyrqm(1) - if Byrqm(2)>Eyrqm(2) - error('Eyrqm(2) must be equal to or greater than Byrqm(2) because of the same year') - end -end - - -Yr = Byrqm(1)+[0:Eyrqm(1)-Byrqm(1)]'; - -if length(Yr)>=2 % there are years and quarters (months) between Byrqm and Eyrqm - n=length(Yr)-2; - C=zeros(n*q_m,2); - C(:,1) = kron(Yr(2:end-1),ones(q_m,1)); - C(:,2) = kron(ones(n,1),[1:q_m]'); - - %* initialize a matrix of years and quarters (months) including Byrqm and Eyrqm - Myrqm = zeros((q_m-Byrqm(2)+1)+Eyrqm(2)+n*q_m,2); - - %* Years in between - n1=q_m-Byrqm(2)+1; - n2=Eyrqm(2); - Myrqm(n1+1:end-n2,:) = C; - %* Beginning year - for k=1:n1 - Myrqm(k,:) = [Byrqm(1) Byrqm(2)+k-1]; - end - %* End year - for k=1:n2 - Myrqm(end-Eyrqm(2)+k,:) = [Eyrqm(1) k]; - end -else %* all the data are in the same calendar year - n1=Eyrqm(2)-Byrqm(2)+1; - Myrqm = zeros(n1,2); - for k=1:n1 - Myrqm(k,:) = [Byrqm(1) Byrqm(2)+k-1]; - end -end - -nMyrqm = size(Myrqm,1); +function [Myrqm,nMyrqm] = fn_calyrqm(q_m,Byrqm,Eyrqm) +% [Myrqm,nMyrqm] = fn_calyrqm(q_m,Byrqm,Eyrqm) +% +% Given the beginning and end years and quarters (months), export a matrix of all years and +% quarters (months) for these years and in between +% +% q_m: 4 if quarterly and 12 if monthly +% Byrqm: [year quarter(month)] -- all integers, the begining year and quarter (month) +% Eyrqm: [year quarter(month)] -- all integers, the end year and quarter (month) +%------------------- +% Myrqm: matrix of all years and quarters (months) between and incl. Byrqm and Eyrqm +% nMyrqm: number of data points incl. Byrqm and Eyrqm +% +% Tao Zha, April 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if ~isempty(find(Byrqm-round(Byrqm))) | (q_m-round(q_m)) | ~isempty(find(Byrqm-round(Byrqm))) + error('argin qm, Byrqm, or Eyrqm must of integer') +elseif Byrqm(1)>Eyrqm(1) + error('Eyrqm(1) must be equal to or greater than Byrqm(1)') +elseif Byrqm(1)==Eyrqm(1) + if Byrqm(2)>Eyrqm(2) + error('Eyrqm(2) must be equal to or greater than Byrqm(2) because of the same year') + end +end + + +Yr = Byrqm(1)+[0:Eyrqm(1)-Byrqm(1)]'; + +if length(Yr)>=2 % there are years and quarters (months) between Byrqm and Eyrqm + n=length(Yr)-2; + C=zeros(n*q_m,2); + C(:,1) = kron(Yr(2:end-1),ones(q_m,1)); + C(:,2) = kron(ones(n,1),[1:q_m]'); + + %* initialize a matrix of years and quarters (months) including Byrqm and Eyrqm + Myrqm = zeros((q_m-Byrqm(2)+1)+Eyrqm(2)+n*q_m,2); + + %* Years in between + n1=q_m-Byrqm(2)+1; + n2=Eyrqm(2); + Myrqm(n1+1:end-n2,:) = C; + %* Beginning year + for k=1:n1 + Myrqm(k,:) = [Byrqm(1) Byrqm(2)+k-1]; + end + %* End year + for k=1:n2 + Myrqm(end-Eyrqm(2)+k,:) = [Eyrqm(1) k]; + end +else %* all the data are in the same calendar year + n1=Eyrqm(2)-Byrqm(2)+1; + Myrqm = zeros(n1,2); + for k=1:n1 + Myrqm(k,:) = [Byrqm(1) Byrqm(2)+k-1]; + end +end + +nMyrqm = size(Myrqm,1); diff --git a/matlab/ms-sbvar/cstz/fn_dataext.m b/matlab/ms-sbvar/cstz/fn_dataext.m index d8d295107f7a836c34722960100a0dc8e7cccf99..7df59f93ef108f62566eed8c3fd4925d50c3b78f 100644 --- a/matlab/ms-sbvar/cstz/fn_dataext.m +++ b/matlab/ms-sbvar/cstz/fn_dataext.m @@ -1,60 +1,60 @@ -function [xdsube,Brow,Erow] = fn_dataext(Byrqm,Eyrqm,xdatae) -% xdsube = dataext(Byrqm,Eyrqm,xdatae) -% Extract subset of xdatae, given Byrqm and Eyrqm -% -% Byrqm: [year quarter(month)]: beginning year and period. If Byqm(2)=0, we get annual data. -% Eyrqm: [year period]: end year and period. If Byrqm(2)=0, it must be Eyrqm(2)=0. -% xdatae: all data (some of which may be NaN) with the first 2 columns indicating years and periods. -%---------- -% xdsube: subset of xdatae. -% Brow: the row number in xdatee that marks the first row of xdsube. -% Erow: the row number in xdatee that marks the last row of xdsube. -% -% Tao Zha, April 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if (Byrqm(2)==0) && (Eyrqm(2)~=0) - error('If annual data, make sure both Byrqm(2) and Eyrqm(2) are zero') -end - -Brow = min(find(xdatae(:,1)==Byrqm(1))); -if isempty(Brow) - error('Byrqm is outside the date range indicated by xdatae(:,1:2)!') -end -if Byrqm(2)>0 - nadt=Byrqm(2)-xdatae(Brow,2); - if nadt<0 - error('Byrqm is outside the date range indicated by xdatae(:,1:2)!') - end - Brow=Brow+nadt; -end -% -Erow = min(find(xdatae(:,1)==Eyrqm(1))); -if isempty(Erow) - error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!') -end -if Eyrqm(2)>0 - nadt2=Eyrqm(2)-xdatae(Erow,2); - if nadt<0 - error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!') - end - Erow=Erow+nadt2; -end -% -xdsube = xdatae(Brow:Erow,:); % with dates +function [xdsube,Brow,Erow] = fn_dataext(Byrqm,Eyrqm,xdatae) +% xdsube = dataext(Byrqm,Eyrqm,xdatae) +% Extract subset of xdatae, given Byrqm and Eyrqm +% +% Byrqm: [year quarter(month)]: beginning year and period. If Byqm(2)=0, we get annual data. +% Eyrqm: [year period]: end year and period. If Byrqm(2)=0, it must be Eyrqm(2)=0. +% xdatae: all data (some of which may be NaN) with the first 2 columns indicating years and periods. +%---------- +% xdsube: subset of xdatae. +% Brow: the row number in xdatee that marks the first row of xdsube. +% Erow: the row number in xdatee that marks the last row of xdsube. +% +% Tao Zha, April 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if (Byrqm(2)==0) && (Eyrqm(2)~=0) + error('If annual data, make sure both Byrqm(2) and Eyrqm(2) are zero') +end + +Brow = min(find(xdatae(:,1)==Byrqm(1))); +if isempty(Brow) + error('Byrqm is outside the date range indicated by xdatae(:,1:2)!') +end +if Byrqm(2)>0 + nadt=Byrqm(2)-xdatae(Brow,2); + if nadt<0 + error('Byrqm is outside the date range indicated by xdatae(:,1:2)!') + end + Brow=Brow+nadt; +end +% +Erow = min(find(xdatae(:,1)==Eyrqm(1))); +if isempty(Erow) + error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!') +end +if Eyrqm(2)>0 + nadt2=Eyrqm(2)-xdatae(Erow,2); + if nadt<0 + error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!') + end + Erow=Erow+nadt2; +end +% +xdsube = xdatae(Brow:Erow,:); % with dates diff --git a/matlab/ms-sbvar/cstz/fn_datana.m b/matlab/ms-sbvar/cstz/fn_datana.m index 331b196e8376abe884354771a475125185320b8b..2da6d59a0ab06373348e46623e537cd1047caadf 100644 --- a/matlab/ms-sbvar/cstz/fn_datana.m +++ b/matlab/ms-sbvar/cstz/fn_datana.m @@ -1,117 +1,117 @@ -function [yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(xdatae,q_m,vlistlog,vlistper,Byrqm,Eyrqm) -% [yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(Byrqm,Eyrqm,xdatae,q_m,vlistlog,vlistper) -% -% Generate prior period, period-to-period-in-last-year, and annual growth rates -% For annual rates, works for both calendar and any annual years, depending on Byrqm and Eyrqm -% If monthly data, we haven't got time to get quarterly growth rates yet. -% -% xdatae: all data (logged levels or interest rates/100, some of which may be NaN) with the first -% 2 columns indicating years and periods. -% q_m: quarter or month period -% vlistlog: sublist for logged variables -% vlistper: sublists for percent variables -% Byrqm: [year quarter(month)]: beginning year and period. If Byqm(2)~=1, we don't get -% calendar annual rates. In other words, the first column of yactyge (which -% indicates years) does not mean calendar years. Byqm(2) must be specified; in other -% words, it must be not set to 0 as in, say, fn_dataext. -% Eyrqm: [year period]: end year and period. Eyqm(2) must be specified; in other words, it -% must be not set to 0 as in, say, fn_dataext. -% NOTE: if no inputs Byrqm and Eyrqm are specified, all growth rates begin at xdatae(1,1:2). -%---------- -% yactyrge: annual growth rates with dates in the first 2 columns. -% yactyre: annual average logged level with dates in the 1st 2 columns. -% yactqmyge: period-to-period-in-last-year annual growth rates with dates in the first 2 columns. -% yactqmge: prior-period annualized growth rates with dates in the first 2 columns. -% yactqme: data (logged levels or interest rates/100) with dates in the first 2 columns. -% Same as xdatae but with Brow:Erow. -% -% Tao Zha, April 2000. - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if size(xdatae,1)<2*q_m - error('We need at least two years of xdatae to get annual rates. Check xdatae!!') -end - -if nargin==4 - Brow=1; Erow=length(xdatae(:,1)); - nyr = floor((Erow-Brow+1)/q_m); - yrsg = [xdatae(q_m+1,1):xdatae(q_m+1,1)+nyr-2]'; % for annual growth later on -else - if Byrqm(2)<1 || Eyrqm(2)<1 - error('This function requires specifying both years and months (quarters) in Byrqm and Eyrqm') - end - - Brow = min(find(xdatae(:,1)==Byrqm(1))); - if isempty(Brow) - error('Byrqm is outside the date range of xdatae(:,1:2)!') - end - nadt=Byrqm(2)-xdatae(Brow,2); - if nadt<0 - error('Byrqm is outside the date range indicated by xdatae(:,1:2)!') - end - Brow=Brow+nadt; - % - Erow = min(find(xdatae(:,1)==Eyrqm(1))); - if isempty(Brow) - error('Eyrqm is outside the date range of xdatae(:,1:2)!') - end - nadt2=Eyrqm(2)-xdatae(Erow,2); - if nadt<0 - error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!') - end - Erow=Erow+nadt2; - - nyr = floor((Erow-Brow+1)/q_m); - yrsg = [Byrqm(1)+1:Byrqm(1)+nyr-1]'; % for annual growth later on, which will - % start at Byrqm(1) instead of Byrqm(1)+1 -end -% -yactqme = xdatae(Brow:Erow,:); % with dates -yactqm = yactqme(:,3:end); % only data - -%======== prior period change (annaluized rate) -yactqmg = yactqm(2:end,:); % start at second period to get growth rate -yactqmg(:,vlistlog) = (yactqm(2:end,vlistlog) - yactqm(1:end-1,vlistlog)) .* q_m; - % monthly, 12*log(1+growth rate), annualized growth rate -yactqmg(:,vlistlog) = 100*(exp(yactqmg(:,vlistlog))-1); -yactqmg(:,vlistper) = 100*yactqmg(:,vlistper); -yactqmge = [yactqme(2:end,1:2) yactqmg]; - -%======== change from the last year -yactqmyg = yactqm(q_m+1:end,:); % start at the last-year period to get growth rate -yactqmyg(:,vlistlog) = (yactqm(q_m+1:end,vlistlog) - yactqm(1:end-q_m,vlistlog)); -yactqmyg(:,vlistlog) = 100*(exp(yactqmyg(:,vlistlog))-1); -yactqmyg(:,vlistper) = 100*yactqmyg(:,vlistper); -yactqmyge = [yactqme(q_m+1:end,1:2) yactqmyg]; - -%======== annual growth rates -nvar = length(xdatae(1,3:end)); -ygmts = yactqm(1:nyr*q_m,:); % converted to the multiplication of q_m -ygmts1 = reshape(ygmts,q_m,nyr,nvar); -ygmts2 = sum(ygmts1,1) ./ q_m; -ygmts3 = reshape(ygmts2,nyr,nvar); % converted to annual average series -% -yactyrg = ygmts3(2:end,:); % start at the last-year period to get growth rate -yactyrg(:,vlistlog) = ygmts3(2:end,vlistlog) - ygmts3(1:end-1,vlistlog); - % annual rate: log(1+growth rate) -yactyrg(:,vlistlog) = 100*(exp(yactyrg(:,vlistlog))-1); -yactyrg(:,vlistper) = 100*yactyrg(:,vlistper); -yactyrge = [yrsg zeros(nyr-1,1) yactyrg]; -yrsg1=[yrsg(1)-1:yrsg(end)]'; -yactyre = [yrsg1 zeros(nyr,1) ygmts3]; +function [yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(xdatae,q_m,vlistlog,vlistper,Byrqm,Eyrqm) +% [yactyrge,yactyre,yactqmyge,yactqmge,yactqme] = fn_datana(Byrqm,Eyrqm,xdatae,q_m,vlistlog,vlistper) +% +% Generate prior period, period-to-period-in-last-year, and annual growth rates +% For annual rates, works for both calendar and any annual years, depending on Byrqm and Eyrqm +% If monthly data, we haven't got time to get quarterly growth rates yet. +% +% xdatae: all data (logged levels or interest rates/100, some of which may be NaN) with the first +% 2 columns indicating years and periods. +% q_m: quarter or month period +% vlistlog: sublist for logged variables +% vlistper: sublists for percent variables +% Byrqm: [year quarter(month)]: beginning year and period. If Byqm(2)~=1, we don't get +% calendar annual rates. In other words, the first column of yactyge (which +% indicates years) does not mean calendar years. Byqm(2) must be specified; in other +% words, it must be not set to 0 as in, say, fn_dataext. +% Eyrqm: [year period]: end year and period. Eyqm(2) must be specified; in other words, it +% must be not set to 0 as in, say, fn_dataext. +% NOTE: if no inputs Byrqm and Eyrqm are specified, all growth rates begin at xdatae(1,1:2). +%---------- +% yactyrge: annual growth rates with dates in the first 2 columns. +% yactyre: annual average logged level with dates in the 1st 2 columns. +% yactqmyge: period-to-period-in-last-year annual growth rates with dates in the first 2 columns. +% yactqmge: prior-period annualized growth rates with dates in the first 2 columns. +% yactqme: data (logged levels or interest rates/100) with dates in the first 2 columns. +% Same as xdatae but with Brow:Erow. +% +% Tao Zha, April 2000. + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if size(xdatae,1)<2*q_m + error('We need at least two years of xdatae to get annual rates. Check xdatae!!') +end + +if nargin==4 + Brow=1; Erow=length(xdatae(:,1)); + nyr = floor((Erow-Brow+1)/q_m); + yrsg = [xdatae(q_m+1,1):xdatae(q_m+1,1)+nyr-2]'; % for annual growth later on +else + if Byrqm(2)<1 || Eyrqm(2)<1 + error('This function requires specifying both years and months (quarters) in Byrqm and Eyrqm') + end + + Brow = min(find(xdatae(:,1)==Byrqm(1))); + if isempty(Brow) + error('Byrqm is outside the date range of xdatae(:,1:2)!') + end + nadt=Byrqm(2)-xdatae(Brow,2); + if nadt<0 + error('Byrqm is outside the date range indicated by xdatae(:,1:2)!') + end + Brow=Brow+nadt; + % + Erow = min(find(xdatae(:,1)==Eyrqm(1))); + if isempty(Brow) + error('Eyrqm is outside the date range of xdatae(:,1:2)!') + end + nadt2=Eyrqm(2)-xdatae(Erow,2); + if nadt<0 + error('Eyrqm is outside the date range indicated by xdatae(:,1:2)!') + end + Erow=Erow+nadt2; + + nyr = floor((Erow-Brow+1)/q_m); + yrsg = [Byrqm(1)+1:Byrqm(1)+nyr-1]'; % for annual growth later on, which will + % start at Byrqm(1) instead of Byrqm(1)+1 +end +% +yactqme = xdatae(Brow:Erow,:); % with dates +yactqm = yactqme(:,3:end); % only data + +%======== prior period change (annaluized rate) +yactqmg = yactqm(2:end,:); % start at second period to get growth rate +yactqmg(:,vlistlog) = (yactqm(2:end,vlistlog) - yactqm(1:end-1,vlistlog)) .* q_m; + % monthly, 12*log(1+growth rate), annualized growth rate +yactqmg(:,vlistlog) = 100*(exp(yactqmg(:,vlistlog))-1); +yactqmg(:,vlistper) = 100*yactqmg(:,vlistper); +yactqmge = [yactqme(2:end,1:2) yactqmg]; + +%======== change from the last year +yactqmyg = yactqm(q_m+1:end,:); % start at the last-year period to get growth rate +yactqmyg(:,vlistlog) = (yactqm(q_m+1:end,vlistlog) - yactqm(1:end-q_m,vlistlog)); +yactqmyg(:,vlistlog) = 100*(exp(yactqmyg(:,vlistlog))-1); +yactqmyg(:,vlistper) = 100*yactqmyg(:,vlistper); +yactqmyge = [yactqme(q_m+1:end,1:2) yactqmyg]; + +%======== annual growth rates +nvar = length(xdatae(1,3:end)); +ygmts = yactqm(1:nyr*q_m,:); % converted to the multiplication of q_m +ygmts1 = reshape(ygmts,q_m,nyr,nvar); +ygmts2 = sum(ygmts1,1) ./ q_m; +ygmts3 = reshape(ygmts2,nyr,nvar); % converted to annual average series +% +yactyrg = ygmts3(2:end,:); % start at the last-year period to get growth rate +yactyrg(:,vlistlog) = ygmts3(2:end,vlistlog) - ygmts3(1:end-1,vlistlog); + % annual rate: log(1+growth rate) +yactyrg(:,vlistlog) = 100*(exp(yactyrg(:,vlistlog))-1); +yactyrg(:,vlistper) = 100*yactyrg(:,vlistper); +yactyrge = [yrsg zeros(nyr-1,1) yactyrg]; +yrsg1=[yrsg(1)-1:yrsg(end)]'; +yactyre = [yrsg1 zeros(nyr,1) ygmts3]; diff --git a/matlab/ms-sbvar/cstz/fn_dataxy.m b/matlab/ms-sbvar/cstz/fn_dataxy.m index 5f9b1bb5f30fe1aaad25384cce95a5ac9d201583..164755312b0462bc1e9d937b6612268d040f187d 100644 --- a/matlab/ms-sbvar/cstz/fn_dataxy.m +++ b/matlab/ms-sbvar/cstz/fn_dataxy.m @@ -1,164 +1,164 @@ -function [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh,e] = fn_dataxy(nvar,lags,z,mu,indxDummy,nexo) -% [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh,e] = fn_dataxy(nvar,lags,z,mu,indxDummy,nexo) -% -% Export arranged data matrices for future estimation of the VAR. -% If indxDummy=0, no mu's are used at all and Bh is the OLS estimate. -% If indxDummy=1, only mu(5) and mu(6) as dummy observations. See fn_rnrprior*.m for using mu(1)-mu(4). -% See Wagonner and Zha's Gibbs sampling paper. -% -% nvar: number of endogenous variables. -% lags: the maximum length of lag -% z: T*(nvar+(nexo-1)) matrix of raw or original data (no manipulation involved) -% with sample size including lags and with exogenous variables other than a constant. -% Order of columns: (1) nvar endogenous variables; (2) (nexo-1) exogenous variables; -% (3) constants will be automatically put in the last column. -% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where -% only mu(5) and mu(6) are used for dummy observations in this function (i.e., -% mu(1)-mu(4) are irrelevant here). See fn_rnrprior*.m for using mu(1)-mu(4). -% mu(1): overall tightness and also for A0; (0.57) -% mu(2): relative tightness for A+; (0.13) -% mu(3): relative tightness for the constant term; (0.1) -% mu(4): tightness on lag decay; (1) -% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5) -% mu(6): weight on single dummy initial observation including constant -% (cointegration, unit roots, and stationarity); (5) -% indxDummy: 1: add dummy observations to the data; 0: no dummy added. -% nexo: number of exogenous variables. The constant term is the default setting. Besides this term, -% we have nexo-1 exogenous variables. Optional. If left blank, nexo is set to 1. -% ------------------- -% xtx: X'X: k-by-k where k=ncoef -% xty: X'Y: k-by-nvar -% yty: Y'Y: nvar-by-nvar -% fss: T: sample size excluding lags. With dummyies, fss=nSample-lags+ndobs. -% phi: X; T-by-k; column: [nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term] -% y: Y: T-by-nvar where T=fss -% ncoef: number of coefficients in *each* equation. RHS coefficients only, nvar*lags+nexo -% xr: the economy size (ncoef-by-ncoef) in qr(phi) so that xr=chol(X'*X) or xr'*xr=X'*X -% Bh: ncoef-by-nvar estimated reduced-form parameter; column: nvar; -% row: ncoef=[nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term] -% e: estimated residual e = y -x*Bh, T-by-nvar -% -% Tao Zha, February 2000. - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if nargin == 5 - nexo=1; % default for constant term -elseif nexo<1 - error('We need at least one exogenous term so nexo must >= 1') -end - -%*** original sample dimension without dummy prior -nSample = size(z,1); % the sample size (including lags, of course) -sb = lags+1; % original beginning without dummies -ncoef = nvar*lags+nexo; % number of coefficients in *each* equation, RHS coefficients only. - -if indxDummy % prior dummy prior - %*** expanded sample dimension by dummy prior - ndobs=nvar+1; % number of dummy observations - fss = nSample+ndobs-lags; - - % - % **** nvar prior dummy observations with the sum of coefficients - % ** construct X for Y = X*B + U where phi = X: (T-lags)*k, Y: (T-lags)*nvar - % ** columns: k = # of [nvar for 1st lag, ..., nvar for last lag, exo var, const] - % ** Now, T=T+ndobs -- added with "ndobs" dummy observations - % - phi = zeros(fss,ncoef); - %* constant term - const = ones(fss,1); - const(1:nvar) = zeros(nvar,1); - phi(:,ncoef) = const; % the first nvar periods: no or zero constant! - %* other exogenous (than) constant term - phi(ndobs+1:end,ncoef-nexo+1:ncoef-1) = z(lags+1:end,nvar+1:nvar+nexo-1); - exox = zeros(ndobs,nexo); - phi(1:ndobs,ncoef-nexo+1:ncoef-1) = exox(:,1:nexo-1); - % this = [] when nexo=1 (no other exogenous than constant) - - xdgel = z(:,1:nvar); % endogenous variable matrix - xdgelint = mean(xdgel(1:lags,:),1); % mean of the first lags initial conditions - %* Dummies - for k=1:nvar - for m=1:lags - phi(ndobs,nvar*(m-1)+k) = xdgelint(k); - phi(k,nvar*(m-1)+k) = xdgelint(k); - % <<>> multiply hyperparameter later - end - end - %* True data - for k=1:lags - phi(ndobs+1:fss,nvar*(k-1)+1:nvar*k) = xdgel(sb-k:nSample-k,:); - % row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, exo var, const] - % Thus, # of columns is nvar*lags+nexo = ncoef. - end - % - % ** Y with "ndobs" dummies added - y = zeros(fss,nvar); - %* Dummies - for k=1:nvar - y(ndobs,k) = xdgelint(k); - y(k,k) = xdgelint(k); - % multiply hyperparameter later - end - %* True data - y(ndobs+1:fss,:) = xdgel(sb:nSample,:); - - phi(1:nvar,:) = 1*mu(5)*phi(1:nvar,:); % standard Sims and Zha prior - y(1:nvar,:) = mu(5)*y(1:nvar,:); % standard Sims and Zha prior - phi(nvar+1,:) = mu(6)*phi(nvar+1,:); - y(nvar+1,:) = mu(6)*y(nvar+1,:); - - [xq,xr]=qr(phi,0); - xtx=xr'*xr; - xty=phi'*y; - [yq,yr]=qr(y,0); - yty=yr'*yr; - Bh = xr\(xr'\xty); % xtx\xty where inv(X'X)*(X'Y) - e=y-phi*Bh; -else - fss = nSample-lags; - % - % ** construct X for Y = X*B + U where phi = X: (T-lags)*k, Y: (T-lags)*nvar - % ** columns: k = # of [nvar for 1st lag, ..., nvar for last lag, exo var, const] - % - phi = zeros(fss,ncoef); - %* constant term - const = ones(fss,1); - phi(:,ncoef) = const; % the first nvar periods: no or zero constant! - %* other exogenous (than) constant term - phi(:,ncoef-nexo+1:ncoef-1) = z(lags+1:end,nvar+1:nvar+nexo-1); - % this = [] when nexo=1 (no other exogenous than constant) - - xdgel = z(:,1:nvar); % endogenous variable matrix - %* True data - for k=1:lags - phi(:,nvar*(k-1)+1:nvar*k) = xdgel(sb-k:nSample-k,:); - % row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, exo var, const] - % Thus, # of columns is nvar*lags+nexo = ncoef. - end - % - y = xdgel(sb:nSample,:); - - [xq,xr]=qr(phi,0); - xtx=xr'*xr; - xty=phi'*y; - [yq,yr]=qr(y,0); - yty=yr'*yr; - Bh = xr\(xr'\xty); % xtx\xty where inv(X'X)*(X'Y) - e=y-phi*Bh; -end +function [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh,e] = fn_dataxy(nvar,lags,z,mu,indxDummy,nexo) +% [xtx,xty,yty,fss,phi,y,ncoef,xr,Bh,e] = fn_dataxy(nvar,lags,z,mu,indxDummy,nexo) +% +% Export arranged data matrices for future estimation of the VAR. +% If indxDummy=0, no mu's are used at all and Bh is the OLS estimate. +% If indxDummy=1, only mu(5) and mu(6) as dummy observations. See fn_rnrprior*.m for using mu(1)-mu(4). +% See Wagonner and Zha's Gibbs sampling paper. +% +% nvar: number of endogenous variables. +% lags: the maximum length of lag +% z: T*(nvar+(nexo-1)) matrix of raw or original data (no manipulation involved) +% with sample size including lags and with exogenous variables other than a constant. +% Order of columns: (1) nvar endogenous variables; (2) (nexo-1) exogenous variables; +% (3) constants will be automatically put in the last column. +% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where +% only mu(5) and mu(6) are used for dummy observations in this function (i.e., +% mu(1)-mu(4) are irrelevant here). See fn_rnrprior*.m for using mu(1)-mu(4). +% mu(1): overall tightness and also for A0; (0.57) +% mu(2): relative tightness for A+; (0.13) +% mu(3): relative tightness for the constant term; (0.1) +% mu(4): tightness on lag decay; (1) +% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5) +% mu(6): weight on single dummy initial observation including constant +% (cointegration, unit roots, and stationarity); (5) +% indxDummy: 1: add dummy observations to the data; 0: no dummy added. +% nexo: number of exogenous variables. The constant term is the default setting. Besides this term, +% we have nexo-1 exogenous variables. Optional. If left blank, nexo is set to 1. +% ------------------- +% xtx: X'X: k-by-k where k=ncoef +% xty: X'Y: k-by-nvar +% yty: Y'Y: nvar-by-nvar +% fss: T: sample size excluding lags. With dummyies, fss=nSample-lags+ndobs. +% phi: X; T-by-k; column: [nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term] +% y: Y: T-by-nvar where T=fss +% ncoef: number of coefficients in *each* equation. RHS coefficients only, nvar*lags+nexo +% xr: the economy size (ncoef-by-ncoef) in qr(phi) so that xr=chol(X'*X) or xr'*xr=X'*X +% Bh: ncoef-by-nvar estimated reduced-form parameter; column: nvar; +% row: ncoef=[nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term] +% e: estimated residual e = y -x*Bh, T-by-nvar +% +% Tao Zha, February 2000. + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if nargin == 5 + nexo=1; % default for constant term +elseif nexo<1 + error('We need at least one exogenous term so nexo must >= 1') +end + +%*** original sample dimension without dummy prior +nSample = size(z,1); % the sample size (including lags, of course) +sb = lags+1; % original beginning without dummies +ncoef = nvar*lags+nexo; % number of coefficients in *each* equation, RHS coefficients only. + +if indxDummy % prior dummy prior + %*** expanded sample dimension by dummy prior + ndobs=nvar+1; % number of dummy observations + fss = nSample+ndobs-lags; + + % + % **** nvar prior dummy observations with the sum of coefficients + % ** construct X for Y = X*B + U where phi = X: (T-lags)*k, Y: (T-lags)*nvar + % ** columns: k = # of [nvar for 1st lag, ..., nvar for last lag, exo var, const] + % ** Now, T=T+ndobs -- added with "ndobs" dummy observations + % + phi = zeros(fss,ncoef); + %* constant term + const = ones(fss,1); + const(1:nvar) = zeros(nvar,1); + phi(:,ncoef) = const; % the first nvar periods: no or zero constant! + %* other exogenous (than) constant term + phi(ndobs+1:end,ncoef-nexo+1:ncoef-1) = z(lags+1:end,nvar+1:nvar+nexo-1); + exox = zeros(ndobs,nexo); + phi(1:ndobs,ncoef-nexo+1:ncoef-1) = exox(:,1:nexo-1); + % this = [] when nexo=1 (no other exogenous than constant) + + xdgel = z(:,1:nvar); % endogenous variable matrix + xdgelint = mean(xdgel(1:lags,:),1); % mean of the first lags initial conditions + %* Dummies + for k=1:nvar + for m=1:lags + phi(ndobs,nvar*(m-1)+k) = xdgelint(k); + phi(k,nvar*(m-1)+k) = xdgelint(k); + % <<>> multiply hyperparameter later + end + end + %* True data + for k=1:lags + phi(ndobs+1:fss,nvar*(k-1)+1:nvar*k) = xdgel(sb-k:nSample-k,:); + % row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, exo var, const] + % Thus, # of columns is nvar*lags+nexo = ncoef. + end + % + % ** Y with "ndobs" dummies added + y = zeros(fss,nvar); + %* Dummies + for k=1:nvar + y(ndobs,k) = xdgelint(k); + y(k,k) = xdgelint(k); + % multiply hyperparameter later + end + %* True data + y(ndobs+1:fss,:) = xdgel(sb:nSample,:); + + phi(1:nvar,:) = 1*mu(5)*phi(1:nvar,:); % standard Sims and Zha prior + y(1:nvar,:) = mu(5)*y(1:nvar,:); % standard Sims and Zha prior + phi(nvar+1,:) = mu(6)*phi(nvar+1,:); + y(nvar+1,:) = mu(6)*y(nvar+1,:); + + [xq,xr]=qr(phi,0); + xtx=xr'*xr; + xty=phi'*y; + [yq,yr]=qr(y,0); + yty=yr'*yr; + Bh = xr\(xr'\xty); % xtx\xty where inv(X'X)*(X'Y) + e=y-phi*Bh; +else + fss = nSample-lags; + % + % ** construct X for Y = X*B + U where phi = X: (T-lags)*k, Y: (T-lags)*nvar + % ** columns: k = # of [nvar for 1st lag, ..., nvar for last lag, exo var, const] + % + phi = zeros(fss,ncoef); + %* constant term + const = ones(fss,1); + phi(:,ncoef) = const; % the first nvar periods: no or zero constant! + %* other exogenous (than) constant term + phi(:,ncoef-nexo+1:ncoef-1) = z(lags+1:end,nvar+1:nvar+nexo-1); + % this = [] when nexo=1 (no other exogenous than constant) + + xdgel = z(:,1:nvar); % endogenous variable matrix + %* True data + for k=1:lags + phi(:,nvar*(k-1)+1:nvar*k) = xdgel(sb-k:nSample-k,:); + % row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, exo var, const] + % Thus, # of columns is nvar*lags+nexo = ncoef. + end + % + y = xdgel(sb:nSample,:); + + [xq,xr]=qr(phi,0); + xtx=xr'*xr; + xty=phi'*y; + [yq,yr]=qr(y,0); + yty=yr'*yr; + Bh = xr\(xr'\xty); % xtx\xty where inv(X'X)*(X'Y) + e=y-phi*Bh; +end diff --git a/matlab/ms-sbvar/cstz/fn_ergodp.m b/matlab/ms-sbvar/cstz/fn_ergodp.m index ef5ac8eb9ff67bd06778ff9e9504c47a84b15e27..7c5c70101f1255838aa2f94227ec76f3423cc59d 100644 --- a/matlab/ms-sbvar/cstz/fn_ergodp.m +++ b/matlab/ms-sbvar/cstz/fn_ergodp.m @@ -1,33 +1,33 @@ -function gpi = fn_ergodp(P) -% gpi = fn_ergodp(P) -% Compute the ergodic probabilities. See Hamilton p.681. -% -% P: n-by-n matrix of transition matrix where all elements in each column sum up to 1. -%----- -% gpi: n-by-1 vector of ergodic probabilities. -% -% Tao Zha August 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -[gpim,gpid] = eig(P); % m: matrix; d: diagonal -[gpidv,gpidvinx] = sort(diag(gpid)); -gpidv = fliplr(gpidv); -gpidvinx = flipud(gpidvinx); -gpim = gpim(:,gpidvinx); -gpi = gpim(:,1)/sum(gpim(:,1)); +function gpi = fn_ergodp(P) +% gpi = fn_ergodp(P) +% Compute the ergodic probabilities. See Hamilton p.681. +% +% P: n-by-n matrix of transition matrix where all elements in each column sum up to 1. +%----- +% gpi: n-by-1 vector of ergodic probabilities. +% +% Tao Zha August 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +[gpim,gpid] = eig(P); % m: matrix; d: diagonal +[gpidv,gpidvinx] = sort(diag(gpid)); +gpidv = fliplr(gpidv); +gpidvinx = flipud(gpidvinx); +gpim = gpim(:,gpidvinx); +gpi = gpim(:,1)/sum(gpim(:,1)); diff --git a/matlab/ms-sbvar/cstz/fn_fcstidcnd.m b/matlab/ms-sbvar/cstz/fn_fcstidcnd.m index 3a15aab0d037866155d7db05361c23400e6f5064..0328793d93fceeb79f89c6a047d0ce96a33ab88e 100644 --- a/matlab/ms-sbvar/cstz/fn_fcstidcnd.m +++ b/matlab/ms-sbvar/cstz/fn_fcstidcnd.m @@ -1,322 +1,322 @@ -function [yhat,Estr,rcon,Rcon,u,v,d] = fn_fcstidcnd(valuecon,stepcon,varcon,nstepsm,... - nconstr,eq_ms,nvar,lags,phil,Aband,Sband,yfore_h,imf3s_h,A0_h,Bh_h,... - forep,TLindx,TLnumber,nCms,eq_Cms) -% [yhat,Estr,rcon,Rcon,u,v,d] = fn_fcstidcnd(valuecon,stepcon,varcon,nstepsm,... -% nconstr,eq_ms,nvar,lags,phil,Aband,Sband,yfore_h,imf3s_h,A0_h,Bh_h,... -% forep,TLindx,TLnumber,nCms,eq_Cms) -% -% Conditional forecasting in the identified model with or without error bands -% It handles conditions on average values as well, so "valuecon" must be -% expressed at average (NOT sum) level. -% Aband is used only once when nconstr>0 and Aband=1, where Gibbs sampler may be used -% Unconditional forecast when imf3s_h, etc is fixed and nconstr=0. -% -% valuecon: vector of values conditioned -% stepcon: sequence (cell) of steps conditioned; if length(stepcon{i}) > 1, the condition -% is then an arithmetic average of log(y) over the stepcon{i} period. -% varcon: vector of variables conditioned -% nconstr: number of DLS constraints -% nstepsm: maximum number of steps in all DLS constraints -% nvar: number of variables in the BVAR model -% lags: number of lags in the BVAR model -% phil: the 1-by-(nvar*lags+1) data matrix where k=nvar*lags+1 -% (last period plus lags before the beginning of forecast) -% Aband: 1: draws from A0 and Bh; 0: no draws -% Sband: 1: draws from random shocks E; 0: no random shocks -% yfore_h: uncondtional forecasts: forep-by-nvar. Never used when nconstr=0. -% In this case, set it to []; -% imf3s_h: 3-dimensional impulse responses matrix: impsteps-by-nvar shocks-by-nvar responses -% Never used when nconstr=0. In this case, set it to []; -% A0_h: A0 contemporaneous parameter matrix -% Bh_h: reduced-form parameter matrix: k-by-nvar, y(t) = X(t)*Bh+e(t) -% where X(t) is k-by-nvar and y(t) is 1-by-nvar -% forep: # of forecast periods (e.g., monthly for a monthly model) -% TLindx: 1-by-nCms vector of 1's and 0's, indicating tight or loose; 1: tighter, 0: looser -% Used only when /* (MS draws) is activated. Right now, MS shocks are deterministic. -% TLnumber: 1-by-nCms vector, lower bound for tight and upper bound for loose -% nCms: # of LZ conditions -% eq_Cms: equation location of MS shocks -% ------ -% yhat: conditional forecasts: forep-by-nvar -% Estr: backed-out structural shocks (from N(0,1)) -% rcon: vector - the difference between valuecon and log(yfore) (unconditional forecasts) -% Rcon: k-by-q (q constranits and k=nvar*max(nsteps)) so that -% Rcon'*e = rcon where e is k-by-1 -% [u,d,v]: svd(Rcon,0) -% -%% See Zha's note "Forecast (1)" p. 5, RATS manual (some errors in RATS), etc. -% -%% Some notations: y(t+1) = y(t)B1 + e(t+1)inv(A0). e(t+1) is 1-by-n. -%% Let r(t+1)=e(t+1)inv(A0) + e(t+2)C + .... where inv(A0) is impulse -%% response at t=1, C at t=2, etc. The row of inv(A0) or C is -%% all responses to one shock. -%% Let r be q-by-1 (such as r(1) = r(t+1) -%% = y(t+1) (constrained) - y(t+1) (forecast)). -%% Use impulse responses to find out R (k-by-q) where k=nvar*nsteps -%% where nsteps the largest constrained step. The key of the program -%% is to creat R using impulse responses -%% Optimal solution for shock e where R'*e=r and e is k-by-1 is -%% e = R*inv(R'*R)*r and k>=q -% -% Copyright (c) March 1998 by Tao Zha. Revised November 1998; -% 3/20/99 Disenabled draws of MS shcoks. To enable it, activate /* part -% 3/20/99 Added A0_h and forep and deleted Cms as input argument. Previous -% programs may not be compatible. -% 3/15/2004 There are some BUG problems when calling fn_fcstcnd.m(). - -% Copyright (C) 1998-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -DLSIdShock = ~isempty(eq_ms); % if not empty, the MS shock is identified as in DLS - -impsteps=size(imf3s_h,1); -if (forep<nstepsm) || (impsteps<nstepsm) - disp('Increase # of forecast or impulse steps!!') - disp('Or decrease # of constraints (nconstr) or constrained steps (stepcon(i))!!') - error('Maximum of conditional steps > # of forecast or impulse steps!!') -end -kts = nvar*nstepsm; % k -- ts: total shocks some of which are restricted and others - % are free. -%*** initializing -Rcon = zeros(kts,nconstr); % R: k-by-q -Econ = zeros(kts,1); % E: k-by-1 -rcon = zeros(nconstr,1); % r: q-by-1 -%rcon=valuecon-diag(yfore(stepcon,varcon)); % another way is to use "loop" below. -tcwc = nvar*lags; % total coefficients without constant -phi=phil; - - - -%---------------------------------------------------- -% Form rcon, Rcon, and Econ (the mean of structural shocks) -%---------------------------------------------------- -if nconstr - A0in = reshape(imf3s_h(1,:,:),nvar,nvar); % nvar shocks-by-nvar responses - for i=1:nconstr - rcon(i)=length(stepcon{i})*valuecon(i) - ... - sum(yfore_h(stepcon{i},varcon(i)),1); %<<>> - Rmat = zeros(nstepsm,nvar); - r2mat = zeros(nstepsm,1); % simply one identified equation - % Must be here inside the loop because it's matrix of one column of Rcon - for j=1:length(stepcon{i}) - if DLSIdShock % Assuming the Fed can't see all other shocks within a month - Rmat(1:stepcon{i}(j),eq_ms) = Rmat(1:stepcon{i}(j),eq_ms) + ... - imf3s_h(stepcon{i}(j):-1:1,eq_ms,varcon(i)); - % Rmat: row--nstepsm, column--nvar shocks (here all shocks except - % the identified one are set to zero) for a particular - % endogenous variable 'varcon(i)'. See Zha Forcast (1), pp.6-7 - else % Rcon random with (A0,A+) - Rmat(1:stepcon{i}(j),:) = Rmat(1:stepcon{i}(j),:) + ... - imf3s_h(stepcon{i}(j):-1:1,:,varcon(i)); - % Rmat: row--nstepsm, column--nvar shocks (here all shocks are - % *not* set to zero) for a particular endogenous - % variable 'varcon(i)'. See Zha Forcast (1), pp.6-7 - end - end - Rmatt = Rmat'; % Now, nvar-by-nstepsm. I think here is where RATS has an error - % i.e. "OVERR" is not transposed when overlaid to "CAPR" - Rcon(:,i)=Rmatt(:); % Rcon: k-by-q where q=nconstr - end - - [u d v]=svd(Rcon,0); %trial - %???? Can we reduce the time by computing inv(R'*R) directly? - % rtr = Rcon'*Rcon; %trial - % rtrinv = inv(Rcon'*Rcon); %trial - vd=v.*(ones(size(v,2),1)*diag(d)'); %trial - dinv = 1./diag(d); % inv(diag(d)) - vdinv=v.*(ones(size(v,2),1)*dinv'); %trial - rtr=vd*vd'; % R'*R - rtrinv = vdinv*vdinv'; % inv(R'*R) - - Econ = Rcon*rtrinv*rcon; % E = R*inv(R'R)*r; the mean of structural shocks -else - Econ = zeros(kts,1); % the mean of shocks is zero under no variable condition - Rcon = NaN; - rcon = NaN; - u = NaN; - d = NaN; - v = NaN; -end - - - -%--------------------------------------- -% No uncertainty at all or only random (A0,A+) -% In other words, no future shocks -%--------------------------------------- -if (~Sband) %|| (nconstr && (length(eq_ms)==1)) - % length(eq_ms)==1 implies one-one mapping between MS shocks and, say, FFR - % if nstepsm==nconstr. If this condition does not hold, this procedure - % is incorrect. I don't have time to fix it now (3/20/99). So I use - % this as a proximation - Estr = reshape(Econ,nvar,nstepsm); - Estr = Estr'; % transpose so that - % Estr: structural shocks. Row--steps, Column--n shocks - Estr = [Estr;zeros(forep-nstepsm,nvar)]; - % Now, forep-by-nvar -- ready for forecasts - Estr(1:nCms,eq_Cms) = TLnumber(:); - Ures = Estr/A0_h; % nstepsm-by-nvar - % Ures: reduced-form residuals. Row--steps; Column--n shocks - - % ** reconstruct x(t) for y(t+h) = x(t+h-1)*B - % ** where phi = x(t+h-1) with last column being constant - % - yhat = zeros(forep,nvar); - for k=1:forep - yhat(k,:) = phi*Bh_h + Ures(k,:); - phi(1,nvar+1:tcwc) = phi(1,1:tcwc-nvar); - phi(1,1:nvar) = yhat(k,:); - % - end - -%--------------------------------------- -% With random future shocks and possibly (A0,A+) depending -% on if imf3s_h is random -%--------------------------------------- -else - %-------------- - % Condition on variables and A random - %-------------- - if nconstr && Aband - warning(' ') - disp('This situation (both E and A random) is still under construction') - disp('It is closely related to Waggoner and Zha ReStat Gibbs sampling method') - disp('Please press ctrl-c to abort') - pause - elseif nconstr - %-------------- - % Condition on variables and DLS MS shock, no A random but S random - %-------------- - if DLSIdShock % other shocks are indepedent of the eq_ms shock - % 3/20/99 The following may be problematic because Osk should depend - % on u (A0_h and Bh_h) in general. I haven't worked out any good version - %/* - % Osk = randn(kts,1); % other shocks - % for j=1:nstepsm - % Osk(nvar*(j-1)+eq_ms)=0; % no shock to the MS or identified equation - % end - % Estr = Econ + Osk; % Econ is non zero only at position - % % eq_ms*j where j=1:nstepsm - % Estr = reshape(Estr,nvar,nstepsm); - % Estr = Estr'; % transpose so that - % % Estr: structural shocks. Row--steps, Column--n shocks - % Estr = [Estr;randn(forep-nstepsm,nvar)]; - % % Now, forep-by-nvar -- ready for forecasts - % - disp('DLS') - Ome = eye(kts) - u*u'; % note, I-u*u' = I - R*inv(R'*R)*R' - %[u1 d1 v1] = svd(Ome); % too slow - [u1 d1] = eig(Ome); - Stdcon = u1*diag(sqrt(diag(abs(d1)))); % lower triagular chol of conditional variance - % see Zha's forecast (1), p.17 - tmp1=zeros(nvar,nstepsm); - tmp1(eq_ms,:)=randn(1,nstepsm); - tmp2=tmp1(:); - %Estr1 = Econ + Stdcon*randn(kts,1); - %jnk = reshape(Stdcon*tmp2,nvar,nstepsm) - Estr1 = Econ + Stdcon*tmp2; - Estr2 = reshape(Estr1,nvar,nstepsm); - Estr2 = Estr2'; % transpose so that - % Estr2: structural shocks. Row--nstepsm, Column--n shocks - Estr = [Estr2;randn(forep-nstepsm,nvar)]; - % Now, forep-by-nvar -- ready for forecasts - else - Ome = eye(kts) - u*u'; % note, I-u*u' = I - R*inv(R'*R)*R' - %[u1 d1 v1] = svd(Ome); % too slow - [u1 d1] = eig(Ome); - Stdcon = u1*diag(sqrt(diag(abs(d1)))); % lower triagular chol of conditional variance - % see Zha's forecast (1), p.17 - %-------------- - % Condition on variables and LZ MS shock, no A random but S random - % This section has not be tested yet, 10/14/98 - %-------------- - if nCms - Estr1 = Econ + Stdcon*randn(kts,1); - Estr2 = reshape(Estr1,nvar,nstepsm); - Estr2 = Estr2'; % transpose so that - % Estr2: structural shocks. Row--nstepsm, Column--n shocks - Estr = [Estr2;randn(forep-nstepsm,nvar)]; - % Now, forep-by-nvar -- ready for forecasts - Estr(1:nCms,eq_Cms) = TLnumber(:); - - %/* draw MS shocks - % for k=1:nCms - % if TLindx(k) % tighter - % while (Estr(k,eq_Cms)<TLnumber(k)) - % Estr(k,eq_Cms) = randn(1,1); - % end - % else % looser - % while (Estr(k,eq_Cms)>TLnumber(k)) - % Estr(k,eq_Cms) = randn(1,1); - % end - % end - % end - %-------------- - % Condition on variables only, no A random but S random - %-------------- - else - Estr1 = Econ + Stdcon*randn(kts,1); - Estr2 = reshape(Estr1,nvar,nstepsm); - Estr2 = Estr2'; % transpose so that - % Estr2: structural shocks. Row--nstepsm, Column--n shocks - Estr = [Estr2;randn(forep-nstepsm,nvar)]; - % Now, forep-by-nvar -- ready for forecasts - end - end - %-------------- - % Condition on LZ MS shocks only, S random and possibly A random depending on - % if A0_h and Bh_h are random - %-------------- - else - if nCms - Estr = randn(forep,nvar); - % Now, forep-by-nvar -- ready for forecasts - Estr(1:nCms,eq_Cms) = TLnumber(:); - - %/* draw MS shocks - % for k=1:nCms - % if TLindx(k) % tighter - % while (Estr(k,eq_Cms)<TLnumber(k)) - % Estr(k,eq_Cms) = randn(1,1); - % end - % else % looser - % while (Estr(k,eq_Cms)>TLnumber(k)) - % Estr(k,eq_Cms) = randn(1,1); - % end - % end - % end - else - Estr = randn(forep,nvar); % Unconditional forecast - % Now, forep-by-nvar -- ready for forecasts - end - end - % - - - Ures = Estr/A0_h; % nstepsm-by-nvar - % Ures: reduced-form residuals. Row--steps; Column--n shocks - - % ** reconstruct x(t) for y(t+h) = x(t+h-1)*B - % ** where phi = x(t+h-1) with last column being constant - % - yhat = zeros(forep,nvar); - for k=1:forep - yhat(k,:) = phi*Bh_h + Ures(k,:); - phi(1,nvar+1:tcwc) = phi(1,1:tcwc-nvar); - phi(1,1:nvar) = yhat(k,:); - end -end +function [yhat,Estr,rcon,Rcon,u,v,d] = fn_fcstidcnd(valuecon,stepcon,varcon,nstepsm,... + nconstr,eq_ms,nvar,lags,phil,Aband,Sband,yfore_h,imf3s_h,A0_h,Bh_h,... + forep,TLindx,TLnumber,nCms,eq_Cms) +% [yhat,Estr,rcon,Rcon,u,v,d] = fn_fcstidcnd(valuecon,stepcon,varcon,nstepsm,... +% nconstr,eq_ms,nvar,lags,phil,Aband,Sband,yfore_h,imf3s_h,A0_h,Bh_h,... +% forep,TLindx,TLnumber,nCms,eq_Cms) +% +% Conditional forecasting in the identified model with or without error bands +% It handles conditions on average values as well, so "valuecon" must be +% expressed at average (NOT sum) level. +% Aband is used only once when nconstr>0 and Aband=1, where Gibbs sampler may be used +% Unconditional forecast when imf3s_h, etc is fixed and nconstr=0. +% +% valuecon: vector of values conditioned +% stepcon: sequence (cell) of steps conditioned; if length(stepcon{i}) > 1, the condition +% is then an arithmetic average of log(y) over the stepcon{i} period. +% varcon: vector of variables conditioned +% nconstr: number of DLS constraints +% nstepsm: maximum number of steps in all DLS constraints +% nvar: number of variables in the BVAR model +% lags: number of lags in the BVAR model +% phil: the 1-by-(nvar*lags+1) data matrix where k=nvar*lags+1 +% (last period plus lags before the beginning of forecast) +% Aband: 1: draws from A0 and Bh; 0: no draws +% Sband: 1: draws from random shocks E; 0: no random shocks +% yfore_h: uncondtional forecasts: forep-by-nvar. Never used when nconstr=0. +% In this case, set it to []; +% imf3s_h: 3-dimensional impulse responses matrix: impsteps-by-nvar shocks-by-nvar responses +% Never used when nconstr=0. In this case, set it to []; +% A0_h: A0 contemporaneous parameter matrix +% Bh_h: reduced-form parameter matrix: k-by-nvar, y(t) = X(t)*Bh+e(t) +% where X(t) is k-by-nvar and y(t) is 1-by-nvar +% forep: # of forecast periods (e.g., monthly for a monthly model) +% TLindx: 1-by-nCms vector of 1's and 0's, indicating tight or loose; 1: tighter, 0: looser +% Used only when /* (MS draws) is activated. Right now, MS shocks are deterministic. +% TLnumber: 1-by-nCms vector, lower bound for tight and upper bound for loose +% nCms: # of LZ conditions +% eq_Cms: equation location of MS shocks +% ------ +% yhat: conditional forecasts: forep-by-nvar +% Estr: backed-out structural shocks (from N(0,1)) +% rcon: vector - the difference between valuecon and log(yfore) (unconditional forecasts) +% Rcon: k-by-q (q constranits and k=nvar*max(nsteps)) so that +% Rcon'*e = rcon where e is k-by-1 +% [u,d,v]: svd(Rcon,0) +% +%% See Zha's note "Forecast (1)" p. 5, RATS manual (some errors in RATS), etc. +% +%% Some notations: y(t+1) = y(t)B1 + e(t+1)inv(A0). e(t+1) is 1-by-n. +%% Let r(t+1)=e(t+1)inv(A0) + e(t+2)C + .... where inv(A0) is impulse +%% response at t=1, C at t=2, etc. The row of inv(A0) or C is +%% all responses to one shock. +%% Let r be q-by-1 (such as r(1) = r(t+1) +%% = y(t+1) (constrained) - y(t+1) (forecast)). +%% Use impulse responses to find out R (k-by-q) where k=nvar*nsteps +%% where nsteps the largest constrained step. The key of the program +%% is to creat R using impulse responses +%% Optimal solution for shock e where R'*e=r and e is k-by-1 is +%% e = R*inv(R'*R)*r and k>=q +% +% Copyright (c) March 1998 by Tao Zha. Revised November 1998; +% 3/20/99 Disenabled draws of MS shcoks. To enable it, activate /* part +% 3/20/99 Added A0_h and forep and deleted Cms as input argument. Previous +% programs may not be compatible. +% 3/15/2004 There are some BUG problems when calling fn_fcstcnd.m(). + +% Copyright (C) 1998-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +DLSIdShock = ~isempty(eq_ms); % if not empty, the MS shock is identified as in DLS + +impsteps=size(imf3s_h,1); +if (forep<nstepsm) || (impsteps<nstepsm) + disp('Increase # of forecast or impulse steps!!') + disp('Or decrease # of constraints (nconstr) or constrained steps (stepcon(i))!!') + error('Maximum of conditional steps > # of forecast or impulse steps!!') +end +kts = nvar*nstepsm; % k -- ts: total shocks some of which are restricted and others + % are free. +%*** initializing +Rcon = zeros(kts,nconstr); % R: k-by-q +Econ = zeros(kts,1); % E: k-by-1 +rcon = zeros(nconstr,1); % r: q-by-1 +%rcon=valuecon-diag(yfore(stepcon,varcon)); % another way is to use "loop" below. +tcwc = nvar*lags; % total coefficients without constant +phi=phil; + + + +%---------------------------------------------------- +% Form rcon, Rcon, and Econ (the mean of structural shocks) +%---------------------------------------------------- +if nconstr + A0in = reshape(imf3s_h(1,:,:),nvar,nvar); % nvar shocks-by-nvar responses + for i=1:nconstr + rcon(i)=length(stepcon{i})*valuecon(i) - ... + sum(yfore_h(stepcon{i},varcon(i)),1); %<<>> + Rmat = zeros(nstepsm,nvar); + r2mat = zeros(nstepsm,1); % simply one identified equation + % Must be here inside the loop because it's matrix of one column of Rcon + for j=1:length(stepcon{i}) + if DLSIdShock % Assuming the Fed can't see all other shocks within a month + Rmat(1:stepcon{i}(j),eq_ms) = Rmat(1:stepcon{i}(j),eq_ms) + ... + imf3s_h(stepcon{i}(j):-1:1,eq_ms,varcon(i)); + % Rmat: row--nstepsm, column--nvar shocks (here all shocks except + % the identified one are set to zero) for a particular + % endogenous variable 'varcon(i)'. See Zha Forcast (1), pp.6-7 + else % Rcon random with (A0,A+) + Rmat(1:stepcon{i}(j),:) = Rmat(1:stepcon{i}(j),:) + ... + imf3s_h(stepcon{i}(j):-1:1,:,varcon(i)); + % Rmat: row--nstepsm, column--nvar shocks (here all shocks are + % *not* set to zero) for a particular endogenous + % variable 'varcon(i)'. See Zha Forcast (1), pp.6-7 + end + end + Rmatt = Rmat'; % Now, nvar-by-nstepsm. I think here is where RATS has an error + % i.e. "OVERR" is not transposed when overlaid to "CAPR" + Rcon(:,i)=Rmatt(:); % Rcon: k-by-q where q=nconstr + end + + [u d v]=svd(Rcon,0); %trial + %???? Can we reduce the time by computing inv(R'*R) directly? + % rtr = Rcon'*Rcon; %trial + % rtrinv = inv(Rcon'*Rcon); %trial + vd=v.*(ones(size(v,2),1)*diag(d)'); %trial + dinv = 1./diag(d); % inv(diag(d)) + vdinv=v.*(ones(size(v,2),1)*dinv'); %trial + rtr=vd*vd'; % R'*R + rtrinv = vdinv*vdinv'; % inv(R'*R) + + Econ = Rcon*rtrinv*rcon; % E = R*inv(R'R)*r; the mean of structural shocks +else + Econ = zeros(kts,1); % the mean of shocks is zero under no variable condition + Rcon = NaN; + rcon = NaN; + u = NaN; + d = NaN; + v = NaN; +end + + + +%--------------------------------------- +% No uncertainty at all or only random (A0,A+) +% In other words, no future shocks +%--------------------------------------- +if (~Sband) %|| (nconstr && (length(eq_ms)==1)) + % length(eq_ms)==1 implies one-one mapping between MS shocks and, say, FFR + % if nstepsm==nconstr. If this condition does not hold, this procedure + % is incorrect. I don't have time to fix it now (3/20/99). So I use + % this as a proximation + Estr = reshape(Econ,nvar,nstepsm); + Estr = Estr'; % transpose so that + % Estr: structural shocks. Row--steps, Column--n shocks + Estr = [Estr;zeros(forep-nstepsm,nvar)]; + % Now, forep-by-nvar -- ready for forecasts + Estr(1:nCms,eq_Cms) = TLnumber(:); + Ures = Estr/A0_h; % nstepsm-by-nvar + % Ures: reduced-form residuals. Row--steps; Column--n shocks + + % ** reconstruct x(t) for y(t+h) = x(t+h-1)*B + % ** where phi = x(t+h-1) with last column being constant + % + yhat = zeros(forep,nvar); + for k=1:forep + yhat(k,:) = phi*Bh_h + Ures(k,:); + phi(1,nvar+1:tcwc) = phi(1,1:tcwc-nvar); + phi(1,1:nvar) = yhat(k,:); + % + end + +%--------------------------------------- +% With random future shocks and possibly (A0,A+) depending +% on if imf3s_h is random +%--------------------------------------- +else + %-------------- + % Condition on variables and A random + %-------------- + if nconstr && Aband + warning(' ') + disp('This situation (both E and A random) is still under construction') + disp('It is closely related to Waggoner and Zha ReStat Gibbs sampling method') + disp('Please press ctrl-c to abort') + pause + elseif nconstr + %-------------- + % Condition on variables and DLS MS shock, no A random but S random + %-------------- + if DLSIdShock % other shocks are indepedent of the eq_ms shock + % 3/20/99 The following may be problematic because Osk should depend + % on u (A0_h and Bh_h) in general. I haven't worked out any good version + %/* + % Osk = randn(kts,1); % other shocks + % for j=1:nstepsm + % Osk(nvar*(j-1)+eq_ms)=0; % no shock to the MS or identified equation + % end + % Estr = Econ + Osk; % Econ is non zero only at position + % % eq_ms*j where j=1:nstepsm + % Estr = reshape(Estr,nvar,nstepsm); + % Estr = Estr'; % transpose so that + % % Estr: structural shocks. Row--steps, Column--n shocks + % Estr = [Estr;randn(forep-nstepsm,nvar)]; + % % Now, forep-by-nvar -- ready for forecasts + % + disp('DLS') + Ome = eye(kts) - u*u'; % note, I-u*u' = I - R*inv(R'*R)*R' + %[u1 d1 v1] = svd(Ome); % too slow + [u1 d1] = eig(Ome); + Stdcon = u1*diag(sqrt(diag(abs(d1)))); % lower triagular chol of conditional variance + % see Zha's forecast (1), p.17 + tmp1=zeros(nvar,nstepsm); + tmp1(eq_ms,:)=randn(1,nstepsm); + tmp2=tmp1(:); + %Estr1 = Econ + Stdcon*randn(kts,1); + %jnk = reshape(Stdcon*tmp2,nvar,nstepsm) + Estr1 = Econ + Stdcon*tmp2; + Estr2 = reshape(Estr1,nvar,nstepsm); + Estr2 = Estr2'; % transpose so that + % Estr2: structural shocks. Row--nstepsm, Column--n shocks + Estr = [Estr2;randn(forep-nstepsm,nvar)]; + % Now, forep-by-nvar -- ready for forecasts + else + Ome = eye(kts) - u*u'; % note, I-u*u' = I - R*inv(R'*R)*R' + %[u1 d1 v1] = svd(Ome); % too slow + [u1 d1] = eig(Ome); + Stdcon = u1*diag(sqrt(diag(abs(d1)))); % lower triagular chol of conditional variance + % see Zha's forecast (1), p.17 + %-------------- + % Condition on variables and LZ MS shock, no A random but S random + % This section has not be tested yet, 10/14/98 + %-------------- + if nCms + Estr1 = Econ + Stdcon*randn(kts,1); + Estr2 = reshape(Estr1,nvar,nstepsm); + Estr2 = Estr2'; % transpose so that + % Estr2: structural shocks. Row--nstepsm, Column--n shocks + Estr = [Estr2;randn(forep-nstepsm,nvar)]; + % Now, forep-by-nvar -- ready for forecasts + Estr(1:nCms,eq_Cms) = TLnumber(:); + + %/* draw MS shocks + % for k=1:nCms + % if TLindx(k) % tighter + % while (Estr(k,eq_Cms)<TLnumber(k)) + % Estr(k,eq_Cms) = randn(1,1); + % end + % else % looser + % while (Estr(k,eq_Cms)>TLnumber(k)) + % Estr(k,eq_Cms) = randn(1,1); + % end + % end + % end + %-------------- + % Condition on variables only, no A random but S random + %-------------- + else + Estr1 = Econ + Stdcon*randn(kts,1); + Estr2 = reshape(Estr1,nvar,nstepsm); + Estr2 = Estr2'; % transpose so that + % Estr2: structural shocks. Row--nstepsm, Column--n shocks + Estr = [Estr2;randn(forep-nstepsm,nvar)]; + % Now, forep-by-nvar -- ready for forecasts + end + end + %-------------- + % Condition on LZ MS shocks only, S random and possibly A random depending on + % if A0_h and Bh_h are random + %-------------- + else + if nCms + Estr = randn(forep,nvar); + % Now, forep-by-nvar -- ready for forecasts + Estr(1:nCms,eq_Cms) = TLnumber(:); + + %/* draw MS shocks + % for k=1:nCms + % if TLindx(k) % tighter + % while (Estr(k,eq_Cms)<TLnumber(k)) + % Estr(k,eq_Cms) = randn(1,1); + % end + % else % looser + % while (Estr(k,eq_Cms)>TLnumber(k)) + % Estr(k,eq_Cms) = randn(1,1); + % end + % end + % end + else + Estr = randn(forep,nvar); % Unconditional forecast + % Now, forep-by-nvar -- ready for forecasts + end + end + % + + + Ures = Estr/A0_h; % nstepsm-by-nvar + % Ures: reduced-form residuals. Row--steps; Column--n shocks + + % ** reconstruct x(t) for y(t+h) = x(t+h-1)*B + % ** where phi = x(t+h-1) with last column being constant + % + yhat = zeros(forep,nvar); + for k=1:forep + yhat(k,:) = phi*Bh_h + Ures(k,:); + phi(1,nvar+1:tcwc) = phi(1,1:tcwc-nvar); + phi(1,1:nvar) = yhat(k,:); + end +end diff --git a/matlab/ms-sbvar/cstz/fn_forecast.m b/matlab/ms-sbvar/cstz/fn_forecast.m index 3efaf0d641558188d978f1b22685ddc50a759135..644131d7574ac22662d8c060f98c40d81138a217 100644 --- a/matlab/ms-sbvar/cstz/fn_forecast.m +++ b/matlab/ms-sbvar/cstz/fn_forecast.m @@ -1,74 +1,74 @@ -function yhat = fn_forecast(Bh,phi,nn,nexo,Xfexo) -% yhat = fn_forecast(Bh,phi,nn,nexo,Xfexo) -% Unconditional forecating without shocks. -% y_hat(t+h) = c + x_hat(t+h-1)*Bh, X: 1*k; Bh: k*nvar; y_hat: 1*nvar -% -% Bh: k-by-nvar, the (posterior) estimate of B. -% phi: the 1-by-(nvar*lags+nexo) data matrix X where k=nvar*lags+1 -% (last period plus lags before the beginning of forecast). -% nn: [nvar,lags,nfqm], nfqm: forecast periods (months or quarters). -% nexo: number of exogenous variables. The constant term is the default setting. -% Besides this term, we have nexo-1 exogenous variables. -% Xfexo: nfqm-by-nexo-1 vector of exoenous variables in the forecast horizon where -% nfqm: number of forecast periods. -%----------- -% yhat: nfqm-by-nvar forecast. -% -% See fn_forecastsim.m with shocks; fn_forecaststre.m. - -% Copyright (C) 2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if nargin == 3 - nexo=1; % default for constant term -elseif nexo<1 - error('We need at least one exogenous term so nexo must >= 1') -end - -% ** setup -nvar = nn(1); -lags = nn(2); -nfqm = nn(3); -tcwx = nvar*lags; % total coefficeint without exogenous variables - -if nexo>1 - if (nfqm > size(Xfexo,1)) - disp(' ') - warning('Make sure the forecast horizon in the exogenous variable matrix Xfexo > forecast periods') - disp('Press ctrl-c to abort') - pause - elseif ((nexo-1) ~= size(Xfexo,2)) - disp(' ') - warning('Make sure that nexo matchs the exogenous variable matrix Xfexo') - disp('Press ctrl-c to abort') - pause - end -end - -% ** reconstruct x(t) for y(t+h) = x(t+h-1)*B -% ** where phi = x(t+h-1) with last column being constant -yhat = zeros(nfqm,nvar); -for k=1:nfqm - yhat(k,:) = phi*Bh; - %*** lagged endogenous variables - phi(1,nvar+1:tcwx) = phi(1,1:tcwx-nvar); - phi(1,1:nvar) = yhat(k,:); - %*** exogenous variables excluding constant terms - if (nexo>1) - phi(1,tcwx+1:tcwx+nexo-1) = Xfexo(k,:); - end -end +function yhat = fn_forecast(Bh,phi,nn,nexo,Xfexo) +% yhat = fn_forecast(Bh,phi,nn,nexo,Xfexo) +% Unconditional forecating without shocks. +% y_hat(t+h) = c + x_hat(t+h-1)*Bh, X: 1*k; Bh: k*nvar; y_hat: 1*nvar +% +% Bh: k-by-nvar, the (posterior) estimate of B. +% phi: the 1-by-(nvar*lags+nexo) data matrix X where k=nvar*lags+1 +% (last period plus lags before the beginning of forecast). +% nn: [nvar,lags,nfqm], nfqm: forecast periods (months or quarters). +% nexo: number of exogenous variables. The constant term is the default setting. +% Besides this term, we have nexo-1 exogenous variables. +% Xfexo: nfqm-by-nexo-1 vector of exoenous variables in the forecast horizon where +% nfqm: number of forecast periods. +%----------- +% yhat: nfqm-by-nvar forecast. +% +% See fn_forecastsim.m with shocks; fn_forecaststre.m. + +% Copyright (C) 2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if nargin == 3 + nexo=1; % default for constant term +elseif nexo<1 + error('We need at least one exogenous term so nexo must >= 1') +end + +% ** setup +nvar = nn(1); +lags = nn(2); +nfqm = nn(3); +tcwx = nvar*lags; % total coefficeint without exogenous variables + +if nexo>1 + if (nfqm > size(Xfexo,1)) + disp(' ') + warning('Make sure the forecast horizon in the exogenous variable matrix Xfexo > forecast periods') + disp('Press ctrl-c to abort') + pause + elseif ((nexo-1) ~= size(Xfexo,2)) + disp(' ') + warning('Make sure that nexo matchs the exogenous variable matrix Xfexo') + disp('Press ctrl-c to abort') + pause + end +end + +% ** reconstruct x(t) for y(t+h) = x(t+h-1)*B +% ** where phi = x(t+h-1) with last column being constant +yhat = zeros(nfqm,nvar); +for k=1:nfqm + yhat(k,:) = phi*Bh; + %*** lagged endogenous variables + phi(1,nvar+1:tcwx) = phi(1,1:tcwx-nvar); + phi(1,1:nvar) = yhat(k,:); + %*** exogenous variables excluding constant terms + if (nexo>1) + phi(1,tcwx+1:tcwx+nexo-1) = Xfexo(k,:); + end +end diff --git a/matlab/ms-sbvar/cstz/fn_foregraph.m b/matlab/ms-sbvar/cstz/fn_foregraph.m index 73596332edb1482c83994b13b38dab6db66b3e20..5aa924fa0504a5e52faa0ad517bb3a95375f0d88 100644 --- a/matlab/ms-sbvar/cstz/fn_foregraph.m +++ b/matlab/ms-sbvar/cstz/fn_foregraph.m @@ -1,66 +1,66 @@ -function fn_foregraph(yfore,yacte,keyindx,rnum,cnum,q_m,ylab,forelabel,conlab) -% -% Graph annual (or calendar year) forecasts vs actual (all data from "msstart.m") -% -% yfore: actual and forecast annual growth data with dates. -% yacte: actual annual growth data with dates. -% keyindx: index for the variables to be graphed -% rnum: number of rows in subplot -% cnum: number of columns in subplot -% q_m: if 4 or 12, quarterly or monthly data -% ylab: string array for the length(keyindx)-by-1 variables -% forelabel: title label for as of time of forecast -% conlab: label for what conditions imposed; e.g., conlab = 'All bar MS shocks inspl' -%------------- -% No output argument for this graph file -% See fn_seriesgraph.m, fn_forerrgraph.m. -% -% Tao Zha, March 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -vyrs = yfore(:,1); -hornum = cell(length(vyrs),1); % horizontal year (number) -count=0; -for k=vyrs' - count=count+1; - jnk=num2str(k); - hornum{count}=jnk(3:4); % e.g., with '1990', we have '90' -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 +function fn_foregraph(yfore,yacte,keyindx,rnum,cnum,q_m,ylab,forelabel,conlab) +% +% Graph annual (or calendar year) forecasts vs actual (all data from "msstart.m") +% +% yfore: actual and forecast annual growth data with dates. +% yacte: actual annual growth data with dates. +% keyindx: index for the variables to be graphed +% rnum: number of rows in subplot +% cnum: number of columns in subplot +% q_m: if 4 or 12, quarterly or monthly data +% ylab: string array for the length(keyindx)-by-1 variables +% forelabel: title label for as of time of forecast +% conlab: label for what conditions imposed; e.g., conlab = 'All bar MS shocks inspl' +%------------- +% No output argument for this graph file +% See fn_seriesgraph.m, fn_forerrgraph.m. +% +% Tao Zha, March 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +vyrs = yfore(:,1); +hornum = cell(length(vyrs),1); % horizontal year (number) +count=0; +for k=vyrs' + count=count+1; + jnk=num2str(k); + hornum{count}=jnk(3:4); % e.g., with '1990', we have '90' +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 index 79fb01524dd5bf12cd354a5196d8ee11d01c03ca..536cac6cc6fd3d0055b6e4d83b0d656454b90fb7 100644 --- a/matlab/ms-sbvar/cstz/fn_fprintmatrix.m +++ b/matlab/ms-sbvar/cstz/fn_fprintmatrix.m @@ -1,60 +1,60 @@ -function fn_fprintmatrix(fid, M, nrows, ncols, indxFloat) -% Prints the matrix to an ascii file indexed by fid. -% -% Inputs: -% fid: Ascii file id. Example: fid = fopen('outdatainp_3s_stv_tvms6lags.prn','a'); -% M: The matrix to be written to the file. -% nrows: Number of rows of M. -% ncols: Number of columns of M. -% indxFloat: 1 if double; -% 2 if single; -% 3 if only 3 significant digits -% 0 if integer. -% - -% Copyright (C) 2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if nrows~=size(M,1) - nrows - size(M,1) - error('fn_fprintmatrix(): Make sure the row number supplied match that of the matrix'); -end -if ncols~=size(M,2) - ncols - size(M,2) - error('fn_fprintmatrix(): Make sure the column number supplied match that of the matrix'); -end -for ki=1:nrows - for kj=1:ncols - if (indxFloat == 1) - fprintf(fid,' %.16e ',M((kj-1)*nrows+ki)); - elseif (indxFloat == 2) - fprintf(fid,' %.8e ',M((kj-1)*nrows+ki)); - elseif (indxFloat == 3) - fprintf(fid,' %.3e ',M((kj-1)*nrows+ki)); - else - fprintf(fid,' %d ',M((kj-1)*nrows+ki)); - end - if (kj==ncols) - fprintf(fid,'\n'); - end - end - if (ki==nrows) - fprintf(fid,'\n\n'); - end -end +function fn_fprintmatrix(fid, M, nrows, ncols, indxFloat) +% Prints the matrix to an ascii file indexed by fid. +% +% Inputs: +% fid: Ascii file id. Example: fid = fopen('outdatainp_3s_stv_tvms6lags.prn','a'); +% M: The matrix to be written to the file. +% nrows: Number of rows of M. +% ncols: Number of columns of M. +% indxFloat: 1 if double; +% 2 if single; +% 3 if only 3 significant digits +% 0 if integer. +% + +% Copyright (C) 2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if nrows~=size(M,1) + nrows + size(M,1) + error('fn_fprintmatrix(): Make sure the row number supplied match that of the matrix'); +end +if ncols~=size(M,2) + ncols + size(M,2) + error('fn_fprintmatrix(): Make sure the column number supplied match that of the matrix'); +end +for ki=1:nrows + for kj=1:ncols + if (indxFloat == 1) + fprintf(fid,' %.16e ',M((kj-1)*nrows+ki)); + elseif (indxFloat == 2) + fprintf(fid,' %.8e ',M((kj-1)*nrows+ki)); + elseif (indxFloat == 3) + fprintf(fid,' %.3e ',M((kj-1)*nrows+ki)); + else + fprintf(fid,' %d ',M((kj-1)*nrows+ki)); + end + if (kj==ncols) + fprintf(fid,'\n'); + end + end + if (ki==nrows) + fprintf(fid,'\n\n'); + end +end diff --git a/matlab/ms-sbvar/cstz/fn_gfmean.m b/matlab/ms-sbvar/cstz/fn_gfmean.m index 7a869f9109fda5e31d0eb290db00e281d2b0edaa..f9c3dd4db0acc8a9a222361b5d531e26152ea1d2 100644 --- a/matlab/ms-sbvar/cstz/fn_gfmean.m +++ b/matlab/ms-sbvar/cstz/fn_gfmean.m @@ -1,58 +1,58 @@ -function [Fmat,gvec] = fn_gfmean(b,P,Vi,nvar,ncoef,n0,np) -% [Fmat,gvec] = fn_gfmean(b,P,Vi,nvar,ncoef,n0,np) -% -% Mean of free lagged parameters g and original lagged parameters F, conditional on comtemporaneous b's -% See Waggoner and Zha's Gibbs sampling -% -% b: sum(n0)-element vector of mean estimate of A0 free parameters -% P: cell(nvar,1). In each cell, the transformation matrix that affects the posterior mean of A+ conditional on A0. -% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith -% equation lagged restriction matrix where k is a total of exogenous variables and -% ri is the number of free parameters. With this transformation, we have fi = Vi*gi -% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a -% vector of free parameters. There must be at least one free parameter left for -% the ith equation. -% nvar: number of endogeous variables -% ncoef: number of original lagged variables per equation -% n0: nvar-element vector, ith element represents the number of free A0 parameters in ith equation -% np: nvar-element vector, ith element represents the number of free A+ parameters in ith equation -%--------------- -% Fmat: ncoef-by-nvar matrix of original lagged parameters A+. Column corresponding to equation. -% gvec: sum(np)-by-1 stacked vector of all free lagged parameters A+. -% -% Tao Zha, February 2000. Revised, August 2000. - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -b=b(:); n0=n0(:); np=np(:); - -n0cum = [0;cumsum(n0)]; -npcum = [0;cumsum(np)]; -gvec = zeros(npcum(end),1); -Fmat = zeros(ncoef,nvar); % ncoef: maximum original lagged parameters per equation - -if ~(length(b)==n0cum(end)) - error('Make inputs n0 and length(b) match exactly') -end - -for kj=1:nvar - bj = b(n0cum(kj)+1:n0cum(kj+1)); - gj = P{kj}*bj; - gvec(npcum(kj)+1:npcum(kj+1)) = gj; - Fmat(:,kj) = Vi{kj}*gj; -end +function [Fmat,gvec] = fn_gfmean(b,P,Vi,nvar,ncoef,n0,np) +% [Fmat,gvec] = fn_gfmean(b,P,Vi,nvar,ncoef,n0,np) +% +% Mean of free lagged parameters g and original lagged parameters F, conditional on comtemporaneous b's +% See Waggoner and Zha's Gibbs sampling +% +% b: sum(n0)-element vector of mean estimate of A0 free parameters +% P: cell(nvar,1). In each cell, the transformation matrix that affects the posterior mean of A+ conditional on A0. +% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith +% equation lagged restriction matrix where k is a total of exogenous variables and +% ri is the number of free parameters. With this transformation, we have fi = Vi*gi +% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a +% vector of free parameters. There must be at least one free parameter left for +% the ith equation. +% nvar: number of endogeous variables +% ncoef: number of original lagged variables per equation +% n0: nvar-element vector, ith element represents the number of free A0 parameters in ith equation +% np: nvar-element vector, ith element represents the number of free A+ parameters in ith equation +%--------------- +% Fmat: ncoef-by-nvar matrix of original lagged parameters A+. Column corresponding to equation. +% gvec: sum(np)-by-1 stacked vector of all free lagged parameters A+. +% +% Tao Zha, February 2000. Revised, August 2000. + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +b=b(:); n0=n0(:); np=np(:); + +n0cum = [0;cumsum(n0)]; +npcum = [0;cumsum(np)]; +gvec = zeros(npcum(end),1); +Fmat = zeros(ncoef,nvar); % ncoef: maximum original lagged parameters per equation + +if ~(length(b)==n0cum(end)) + error('Make inputs n0 and length(b) match exactly') +end + +for kj=1:nvar + bj = b(n0cum(kj)+1:n0cum(kj+1)); + gj = P{kj}*bj; + gvec(npcum(kj)+1:npcum(kj+1)) = gj; + Fmat(:,kj) = Vi{kj}*gj; +end diff --git a/matlab/ms-sbvar/cstz/fn_gibbsrvar.m b/matlab/ms-sbvar/cstz/fn_gibbsrvar.m index e19329062b0d723e9e3edc6bbf2ca9b71b029983..cc5cd075c3996783e012017dbfb0ab5ad68210a3 100644 --- a/matlab/ms-sbvar/cstz/fn_gibbsrvar.m +++ b/matlab/ms-sbvar/cstz/fn_gibbsrvar.m @@ -1,83 +1,83 @@ -function [A0gbs, Wcell] = fn_gibbsrvar(A0gbs,UT,nvar,fss,n0,Indxcol) -% [A0gbs, Wcell] = fn_gibbsrvar(A0gbs,UT,nvar,fss,n0,Indxcol) -% One-step Gibbs sampler for restricted VARs in the structural form (including over-identified cases). -% Reference: "A Gibbs sampler for structural VARs" by D.F. Waggoner and T. Zha, `` -% Journal of Economic Dynamics & Control (JEDC) 28 (2003) 349-366. -% See Note Forecast (2) pp. 44-51, 70-71, and Theorem 1 and Section 3.1 in the WZ JEDC paper. -% -% A0gbs: the last draw of A0 matrix -% UT: cell(nvar,1) -- U_i*T_i in the proof of Theorem 1 where -% (1) a_i = U_i*b_i with b_i being a vector of free parameters -% (2) T_i (q_i-by-q_i) is from T_i*T_i'= S_i = inv(H0inv{i}/T) where H0inv is the inverse of -% the covariance martrix NOT divided by fss and S_i is defined in (14) on p.355 of the WZ JEDC paper. -% nvar: rank of A0 or # of variables -% fss: effective sample size == nSample (T)-lags+# of dummy observations -% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation -% Indxcol: a row vector indicating random columns this Gibbs draws. -% When this input is not supplied, the Gibbs draws all columns -%------------------ -% A0bgs: nvar-by-nvar. New draw of A0 matrix in this Gibbs step -% Wcell: cell(nvar,1). In each cell, columns of Wcell{i} form an orthonormal basis w_1,...,w_qi in Section 3.1 in the WZ paper. Added 9/04. -% -% Written by Tao Zha, August 2000. Revised, September 2004. -% Copyright (c) by Waggoner and Zha - -% Copyright (C) 2000-2011 Tao Zha and Daniel Waggoner -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if (nargin==5), Indxcol=[1:nvar]; end - -%---------------- Local loop for Gibbs given last A0gbs ---------- -Wcell = cell(length(Indxcol),1); -w = zeros(nvar,1); -for ki=Indxcol % given last A0gbs and generate new A0bgs - X = A0gbs; % WZ's Section 4.3 - X(:,ki) = 0; % want to find non-zero sw s.t., X'*w=0 - - - %*** Solving for w and getting an orthonormal basis. See Theorem 1 and also p.48 in Forecast II - [jL,Ux] = lu(X'); - jIx0 = min(find(abs(diag(Ux))<eps)); % if isempty(jIx0), then something is wrong here - w(jIx0+1:end) = 0; % if jIx0+1>end, no effect on w0 - w(jIx0) = 1; - jA = Ux(1:jIx0-1,1:jIx0-1); - jb = Ux(1:jIx0-1,jIx0); - jy = -jA\jb; - w(1:jIx0-1) = jy; - % Note: if jIx0=1 (which almost never happens for numerical stability reasons), no effect on w. - - %*** Constructing orthonormal basis w_1, ..., w_qi at each Gibbs step - w0 = UT{ki}'*w; - w1 = w0/sqrt(sum(w0.^2)); - [W,jnk] = qr(w1); % Columns of W form an orthonormal basis w_1,...,w_qi in Section 3.1 in the WZ paper - - %*** Draw beta's according to Theorem 1 in the WZ JEDC paper. - gkbeta = zeros(n0(ki),1); % qi-by-1: greak beta's - jstd = sqrt(1/fss); - gkbeta(2:end) = jstd*randn(n0(ki)-1,1); % for beta_2, ..., beta_qi - %--- Unnormalized (i.e., not normalized) gamma or 1-d Wishart draw of beta_1 in Theorem 1 of the WZ JEDC paper. - jr = jstd*randn(fss+1,1); - if rand(1)<0.5 - gkbeta(1) = sqrt(jr'*jr); - else - gkbeta(1) = -sqrt(jr'*jr); - end - - %*** Getting a new ki_th column in A0 - A0gbs(:,ki) = UT{ki}*(W*gkbeta); % n-by-1: U_i*T_i*W*beta; - Wcell{ki} = W; %q_i-by-1. -end +function [A0gbs, Wcell] = fn_gibbsrvar(A0gbs,UT,nvar,fss,n0,Indxcol) +% [A0gbs, Wcell] = fn_gibbsrvar(A0gbs,UT,nvar,fss,n0,Indxcol) +% One-step Gibbs sampler for restricted VARs in the structural form (including over-identified cases). +% Reference: "A Gibbs sampler for structural VARs" by D.F. Waggoner and T. Zha, `` +% Journal of Economic Dynamics & Control (JEDC) 28 (2003) 349-366. +% See Note Forecast (2) pp. 44-51, 70-71, and Theorem 1 and Section 3.1 in the WZ JEDC paper. +% +% A0gbs: the last draw of A0 matrix +% UT: cell(nvar,1) -- U_i*T_i in the proof of Theorem 1 where +% (1) a_i = U_i*b_i with b_i being a vector of free parameters +% (2) T_i (q_i-by-q_i) is from T_i*T_i'= S_i = inv(H0inv{i}/T) where H0inv is the inverse of +% the covariance martrix NOT divided by fss and S_i is defined in (14) on p.355 of the WZ JEDC paper. +% nvar: rank of A0 or # of variables +% fss: effective sample size == nSample (T)-lags+# of dummy observations +% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation +% Indxcol: a row vector indicating random columns this Gibbs draws. +% When this input is not supplied, the Gibbs draws all columns +%------------------ +% A0bgs: nvar-by-nvar. New draw of A0 matrix in this Gibbs step +% Wcell: cell(nvar,1). In each cell, columns of Wcell{i} form an orthonormal basis w_1,...,w_qi in Section 3.1 in the WZ paper. Added 9/04. +% +% Written by Tao Zha, August 2000. Revised, September 2004. +% Copyright (c) by Waggoner and Zha + +% Copyright (C) 2000-2011 Tao Zha and Daniel Waggoner +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if (nargin==5), Indxcol=[1:nvar]; end + +%---------------- Local loop for Gibbs given last A0gbs ---------- +Wcell = cell(length(Indxcol),1); +w = zeros(nvar,1); +for ki=Indxcol % given last A0gbs and generate new A0bgs + X = A0gbs; % WZ's Section 4.3 + X(:,ki) = 0; % want to find non-zero sw s.t., X'*w=0 + + + %*** Solving for w and getting an orthonormal basis. See Theorem 1 and also p.48 in Forecast II + [jL,Ux] = lu(X'); + jIx0 = min(find(abs(diag(Ux))<eps)); % if isempty(jIx0), then something is wrong here + w(jIx0+1:end) = 0; % if jIx0+1>end, no effect on w0 + w(jIx0) = 1; + jA = Ux(1:jIx0-1,1:jIx0-1); + jb = Ux(1:jIx0-1,jIx0); + jy = -jA\jb; + w(1:jIx0-1) = jy; + % Note: if jIx0=1 (which almost never happens for numerical stability reasons), no effect on w. + + %*** Constructing orthonormal basis w_1, ..., w_qi at each Gibbs step + w0 = UT{ki}'*w; + w1 = w0/sqrt(sum(w0.^2)); + [W,jnk] = qr(w1); % Columns of W form an orthonormal basis w_1,...,w_qi in Section 3.1 in the WZ paper + + %*** Draw beta's according to Theorem 1 in the WZ JEDC paper. + gkbeta = zeros(n0(ki),1); % qi-by-1: greak beta's + jstd = sqrt(1/fss); + gkbeta(2:end) = jstd*randn(n0(ki)-1,1); % for beta_2, ..., beta_qi + %--- Unnormalized (i.e., not normalized) gamma or 1-d Wishart draw of beta_1 in Theorem 1 of the WZ JEDC paper. + jr = jstd*randn(fss+1,1); + if rand(1)<0.5 + gkbeta(1) = sqrt(jr'*jr); + else + gkbeta(1) = -sqrt(jr'*jr); + end + + %*** Getting a new ki_th column in A0 + A0gbs(:,ki) = UT{ki}*(W*gkbeta); % n-by-1: U_i*T_i*W*beta; + Wcell{ki} = W; %q_i-by-1. +end diff --git a/matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m b/matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m index 3f2d86b548da5cf769b96959d5aec51a22b51b5d..54810c4a14cb8831a12ad20aafc45160b734046d 100644 --- a/matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m +++ b/matlab/ms-sbvar/cstz/fn_gibbsrvar_setup.m @@ -1,74 +1,74 @@ -function [Tinv,UT,VHphalf,PU,VPU] = fn_gibbsrvar_setup(H0inv, Ui, Hpinv, Pmat, Vi, nvar, fss) -% [Tinv,UT,VHphalf,PU,VPU] = fn_gibbsrvar_setup.m(H0inv, Ui, Hpinv, Pmat, Vi, fss, nvar) -% Global setup outside the Gibbs loop to be used by fn_gibbsvar(). -% Reference: "A Gibbs sampler for structural VARs" by D.F. Waggoner and T. Zha, `` -% Journal of Economic Dynamics & Control (JEDC) 28 (2003) 349-366. -% See Note Forecast (2) pp. 44-51, 70-71, and Theorem 1 and Section 3.1 in the WZ JEDC paper. -% -% H0inv: cell(nvar,1). Not divided by T yet. In each cell, inverse of posterior covariance matrix H0. -% The exponential term is b_i'*inv(H0)*b_i for the ith equation where b_i = U_i*a0_i. -% It resembles old SpH or Sbd in the exponent term in posterior of A0, but not divided by T yet. -% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith -% equation contemporaneous restriction matrix where qi is the number of free parameters. -% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector -% of total original parameters and bi is a vector of free parameters. When no -% restrictions are imposed, we have Ui = I. There must be at least one free -% parameter left for the ith equation. Imported from dnrprior.m. -% Hpinv: cell(nvar,1). In each cell, posterior inverse of covariance matrix Hp (A+) for the free parameters -% g_i = V_i*A+(:,i) in the ith equation. -% Pmat: cell(nvar,1). In each cell, the transformation matrix that affects the posterior mean of A+ conditional on A0. -% In other words, the posterior mean (of g_i) = Pmat{i}*b_i where g_i is a column vector of free parameters -% of A+(:,i)) given b_i (b_i is a column vector of free parameters of A0(:,i)). -% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith -% equation lagged restriction matrix where k (ncoef) is a total number of RHS variables and -% ri is the number of free parameters. With this transformation, we have fi = Vi*gi -% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a -% vector of free parameters. There must be at least one free parameter left for -% the ith equation. Imported from dnrprior.m. -% nvar: number of endogenous variables or rank of A0. -% fss: effective sample size (in the exponential term) = nSample - lags + ndobs (ndobs = # of dummy observations -% is set to 0 when fn_rnrprior_covres_dobs() is used where dummy observations are included as part of the explicit prior. -%------------- -% Tinv: cell(nvar,1). In each cell, inv(T_i) for T_iT_i'=S_i where S_i is defined on p.355 of the WZ JEDC paper. -% UT: cell(nvar,1). In each cell, U_i*T_i. -% VHphalf: cell(nvar,1). In each cell, V_i*sqrt(Hp_i). -% PU: cell(nvar,1). In each cell, Pmat{i}*U_i where Pmat{i} = P_i defined in (13) on p.353 of the WZ JEDC paper. -% VPU: cell(nvar,1). In each cell, V_i*P_i*U_i -% -% Written by Tao Zha, September 2004. -% Copyright (c) 2004 by Waggoner and Zha - -% Copyright (C) 2004-2011 Tao Zha and Daniel Waggoner -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -%--- For A0. -Tinv = cell(nvar,1); % in each cell, inv(T_i) for T_iT_i'=S_i where S_i is defined on p.355 of the WZ JEDC paper. -UT = cell(nvar,1); % in each cell, U_i*T_i. -%--- For A+. -VHphalf = cell(nvar,1); % in each cell, V_i*sqrt(Hp_i). -PU = cell(nvar,1); % in each cell, Pmat{i}*U_i where Pmat{i} = P_i defined in (13) on p.353 of the WZ JEDC paper. -VPU = cell(nvar,1); % in each cell, V_i*P_i*U_i -% -for ki=1:nvar - %--- For A0. - Tinv{ki} = chol(H0inv{ki}/fss); % Tinv_i'*Tinv_i = inv(S_i) ==> T_i*T_i' = S_i where S_i = H0inv{i}/fss is defined on p.355 of the WZ JEDC paper. - UT{ki} = Ui{ki}/Tinv{ki}; % n-by-qi: U_i*T_i in (14) on p. 255 of the WZ JEDC paper. - %--- For A+. - VHphalf{ki} = Vi{ki}/chol(Hpinv{ki}); % where chol(Hpinv_i)*chol(Hpinv_i)'=Hpinv_i. - PU{ki} = Pmat{ki}*Ui{ki}'; - VPU{ki} = Vi{ki}*PU{ki}; -end +function [Tinv,UT,VHphalf,PU,VPU] = fn_gibbsrvar_setup(H0inv, Ui, Hpinv, Pmat, Vi, nvar, fss) +% [Tinv,UT,VHphalf,PU,VPU] = fn_gibbsrvar_setup.m(H0inv, Ui, Hpinv, Pmat, Vi, fss, nvar) +% Global setup outside the Gibbs loop to be used by fn_gibbsvar(). +% Reference: "A Gibbs sampler for structural VARs" by D.F. Waggoner and T. Zha, `` +% Journal of Economic Dynamics & Control (JEDC) 28 (2003) 349-366. +% See Note Forecast (2) pp. 44-51, 70-71, and Theorem 1 and Section 3.1 in the WZ JEDC paper. +% +% H0inv: cell(nvar,1). Not divided by T yet. In each cell, inverse of posterior covariance matrix H0. +% The exponential term is b_i'*inv(H0)*b_i for the ith equation where b_i = U_i*a0_i. +% It resembles old SpH or Sbd in the exponent term in posterior of A0, but not divided by T yet. +% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith +% equation contemporaneous restriction matrix where qi is the number of free parameters. +% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector +% of total original parameters and bi is a vector of free parameters. When no +% restrictions are imposed, we have Ui = I. There must be at least one free +% parameter left for the ith equation. Imported from dnrprior.m. +% Hpinv: cell(nvar,1). In each cell, posterior inverse of covariance matrix Hp (A+) for the free parameters +% g_i = V_i*A+(:,i) in the ith equation. +% Pmat: cell(nvar,1). In each cell, the transformation matrix that affects the posterior mean of A+ conditional on A0. +% In other words, the posterior mean (of g_i) = Pmat{i}*b_i where g_i is a column vector of free parameters +% of A+(:,i)) given b_i (b_i is a column vector of free parameters of A0(:,i)). +% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith +% equation lagged restriction matrix where k (ncoef) is a total number of RHS variables and +% ri is the number of free parameters. With this transformation, we have fi = Vi*gi +% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a +% vector of free parameters. There must be at least one free parameter left for +% the ith equation. Imported from dnrprior.m. +% nvar: number of endogenous variables or rank of A0. +% fss: effective sample size (in the exponential term) = nSample - lags + ndobs (ndobs = # of dummy observations +% is set to 0 when fn_rnrprior_covres_dobs() is used where dummy observations are included as part of the explicit prior. +%------------- +% Tinv: cell(nvar,1). In each cell, inv(T_i) for T_iT_i'=S_i where S_i is defined on p.355 of the WZ JEDC paper. +% UT: cell(nvar,1). In each cell, U_i*T_i. +% VHphalf: cell(nvar,1). In each cell, V_i*sqrt(Hp_i). +% PU: cell(nvar,1). In each cell, Pmat{i}*U_i where Pmat{i} = P_i defined in (13) on p.353 of the WZ JEDC paper. +% VPU: cell(nvar,1). In each cell, V_i*P_i*U_i +% +% Written by Tao Zha, September 2004. +% Copyright (c) 2004 by Waggoner and Zha + +% Copyright (C) 2004-2011 Tao Zha and Daniel Waggoner +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +%--- For A0. +Tinv = cell(nvar,1); % in each cell, inv(T_i) for T_iT_i'=S_i where S_i is defined on p.355 of the WZ JEDC paper. +UT = cell(nvar,1); % in each cell, U_i*T_i. +%--- For A+. +VHphalf = cell(nvar,1); % in each cell, V_i*sqrt(Hp_i). +PU = cell(nvar,1); % in each cell, Pmat{i}*U_i where Pmat{i} = P_i defined in (13) on p.353 of the WZ JEDC paper. +VPU = cell(nvar,1); % in each cell, V_i*P_i*U_i +% +for ki=1:nvar + %--- For A0. + Tinv{ki} = chol(H0inv{ki}/fss); % Tinv_i'*Tinv_i = inv(S_i) ==> T_i*T_i' = S_i where S_i = H0inv{i}/fss is defined on p.355 of the WZ JEDC paper. + UT{ki} = Ui{ki}/Tinv{ki}; % n-by-qi: U_i*T_i in (14) on p. 255 of the WZ JEDC paper. + %--- For A+. + VHphalf{ki} = Vi{ki}/chol(Hpinv{ki}); % where chol(Hpinv_i)*chol(Hpinv_i)'=Hpinv_i. + PU{ki} = Pmat{ki}*Ui{ki}'; + VPU{ki} = Vi{ki}*PU{ki}; +end diff --git a/matlab/ms-sbvar/cstz/fn_imcgraph.m b/matlab/ms-sbvar/cstz/fn_imcgraph.m index f4f24c3dbbf7864c142b2400291ce5a08d903b36..4eeebaaace83fd696fd7577099ba32f456aaccb6 100644 --- a/matlab/ms-sbvar/cstz/fn_imcgraph.m +++ b/matlab/ms-sbvar/cstz/fn_imcgraph.m @@ -1,124 +1,124 @@ -function scaleout = fn_imcgraph(imf,nvar,imstp,xlab,ylab,indxGimfml,xTick) -% scaleout = fn_imcgraph(imf,nvar,imstp,xlab,ylab,indxGimfml,xTick) -% imcgraph: impulse, c (column: shock 1 to N), graph the ML point impulse response -% -% imf: imstp-by-nvar^2 matrix of impulse responses, column (responses to 1st shock, responses to 2nd shock -% etc), row (impusle steps), -% nvar: number of variables -% imstp: number of steps of impulse responses -% xlab,ylab: labels -% indxGimfml: 1, graph; 0, no graph -% xTick: optional. Eg: [12 24 36]. -%--------------- -% scaleout: column 1 represents maximums; column 2 minimums. Rows: nvar variables. -% -% NOTE: I added "indxGimfml" so this function may not be compatible with programs -% older than 03/06/99, TZ -% -% See imrgraph, fn_imcerrgraph, fn_imc2errgraph, imrerrgraph, fn_gyrfore in RVARcode - -% Copyright (C) 1999-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if nargin < 7, xTick = []; end - -t = 1:imstp; -temp1=zeros(nvar,1); -temp2=zeros(nvar,1); -maxval=zeros(nvar,1); -minval=zeros(nvar,1); -for i = 1:nvar - for j = 1:nvar - temp1(j)=max(imf(:,(j-1)*nvar+i)); - temp2(j)=min(imf(:,(j-1)*nvar+i)); - end - maxval(i)=max(temp1); - minval(i)=min(temp2); -end - -scaleout = [maxval(:) minval(:)]; - -%-------------- -% Column j: Shock 1 to N; Row i: Responses to -%------------- -if indxGimfml - %figure - rowlabel = 1; - for i = 1:nvar % Responses of - columnlabel = 1; - - if minval(i)<0 - if maxval(i)<=0 - yt=[minval(i) 0]; - else - yt=[minval(i) 0 maxval(i)]; - end - else % (minval(i) >=0) - if maxval(i) > 0 - yt=[0 maxval(i)]; - else % (identically zero responses) - yt=[-1 0 1]; - end - end - - - scale=[1 imstp minval(i) maxval(i)]; - for j = 1:nvar % To shocks - k1=(i-1)*nvar+j; - k2=(j-1)*nvar+i; - subplot(nvar,nvar,k1) - plot(t,imf(:,k2),t,zeros(length(imf(:,k2)),1),'r:'); - - set(gca,'XTick',xTick) - set(gca,'YTick',yt) - grid - - axis(scale); % put limits on both axes. - % - % if maxval(i)>minval(i) - % set(gca,'YLim',[minval(i) maxval(i)]) - % end - - - if isempty(xTick) %1 % No numbers on both axes - set(gca,'XTickLabel',' '); - set(gca,'YTickLabel',' '); - else % Put numbers on both axes - if i<nvar - set(gca,'XTickLabelMode','manual','XTickLabel',[]) - end - if j>1 - set(gca,'YTickLabel',' '); - end - end - - - if rowlabel == 1 - %title(['x' num2str(j)]) - %title(eval(['x' num2str(j)])) - title(char(xlab(j))) - end - if columnlabel == 1 - %ylabel(['x' num2str(i)]) - %ylabel(eval(['x' num2str(i)])) - ylabel(char(ylab(i))) - end - columnlabel = 0; - end - rowlabel = 0; - end -end +function scaleout = fn_imcgraph(imf,nvar,imstp,xlab,ylab,indxGimfml,xTick) +% scaleout = fn_imcgraph(imf,nvar,imstp,xlab,ylab,indxGimfml,xTick) +% imcgraph: impulse, c (column: shock 1 to N), graph the ML point impulse response +% +% imf: imstp-by-nvar^2 matrix of impulse responses, column (responses to 1st shock, responses to 2nd shock +% etc), row (impusle steps), +% nvar: number of variables +% imstp: number of steps of impulse responses +% xlab,ylab: labels +% indxGimfml: 1, graph; 0, no graph +% xTick: optional. Eg: [12 24 36]. +%--------------- +% scaleout: column 1 represents maximums; column 2 minimums. Rows: nvar variables. +% +% NOTE: I added "indxGimfml" so this function may not be compatible with programs +% older than 03/06/99, TZ +% +% See imrgraph, fn_imcerrgraph, fn_imc2errgraph, imrerrgraph, fn_gyrfore in RVARcode + +% Copyright (C) 1999-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if nargin < 7, xTick = []; end + +t = 1:imstp; +temp1=zeros(nvar,1); +temp2=zeros(nvar,1); +maxval=zeros(nvar,1); +minval=zeros(nvar,1); +for i = 1:nvar + for j = 1:nvar + temp1(j)=max(imf(:,(j-1)*nvar+i)); + temp2(j)=min(imf(:,(j-1)*nvar+i)); + end + maxval(i)=max(temp1); + minval(i)=min(temp2); +end + +scaleout = [maxval(:) minval(:)]; + +%-------------- +% Column j: Shock 1 to N; Row i: Responses to +%------------- +if indxGimfml + %figure + rowlabel = 1; + for i = 1:nvar % Responses of + columnlabel = 1; + + if minval(i)<0 + if maxval(i)<=0 + yt=[minval(i) 0]; + else + yt=[minval(i) 0 maxval(i)]; + end + else % (minval(i) >=0) + if maxval(i) > 0 + yt=[0 maxval(i)]; + else % (identically zero responses) + yt=[-1 0 1]; + end + end + + + scale=[1 imstp minval(i) maxval(i)]; + for j = 1:nvar % To shocks + k1=(i-1)*nvar+j; + k2=(j-1)*nvar+i; + subplot(nvar,nvar,k1) + plot(t,imf(:,k2),t,zeros(length(imf(:,k2)),1),'r:'); + + set(gca,'XTick',xTick) + set(gca,'YTick',yt) + grid + + axis(scale); % put limits on both axes. + % + % if maxval(i)>minval(i) + % set(gca,'YLim',[minval(i) maxval(i)]) + % end + + + if isempty(xTick) %1 % No numbers on both axes + set(gca,'XTickLabel',' '); + set(gca,'YTickLabel',' '); + else % Put numbers on both axes + if i<nvar + set(gca,'XTickLabelMode','manual','XTickLabel',[]) + end + if j>1 + set(gca,'YTickLabel',' '); + end + end + + + if rowlabel == 1 + %title(['x' num2str(j)]) + %title(eval(['x' num2str(j)])) + title(char(xlab(j))) + end + if columnlabel == 1 + %ylabel(['x' num2str(i)]) + %ylabel(eval(['x' num2str(i)])) + ylabel(char(ylab(i))) + end + columnlabel = 0; + end + rowlabel = 0; + end +end diff --git a/matlab/ms-sbvar/cstz/fn_impulse.m b/matlab/ms-sbvar/cstz/fn_impulse.m index 799dee81cdbcbe97ebc6d35fb1d58ba4ef8bcb32..fb91120e97319d8c7b99d677c83eb6e5b06a3b78 100644 --- a/matlab/ms-sbvar/cstz/fn_impulse.m +++ b/matlab/ms-sbvar/cstz/fn_impulse.m @@ -1,76 +1,76 @@ -function imf = fn_impulse(Bh,swish,nn) -% Computing impulse functions with -% imf = fn_impulse(Bh,swish,nn) -% imf is in a format that is the SAME as in RATS. -% Column: nvar responses to 1st shock, -% nvar responses to 2nd shock, and so on. -% Row: steps of impulse responses. -%----------------- -% Bh is the estimated reduced form coefficient in the form -% Y(T*nvar) = XB + U, X: T*k (may include all exogenous terms), B: k*nvar. -% The matrix form and dimension are the same as "Bh" from the function "sye.m"; -% Column: 1st lag (with nvar variables) to lags (with nvar variables) + const = k. -% Note: columns correspond to equations. -% swish is the inv(A0) in the structural model y(t)A0 = e(t). -% Note: columns corresponding to equations. -% nn is the numbers of inputs [nvar,lags,# of steps of impulse responses]. -% -% Written by Tao Zha. -% Copyright (c) 1994 by Tao Zha - -% Copyright (C) 1994-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -nvar = nn(1); -lags = nn(2); -imstep = nn(3); % number of steps for impulse responses - -Ah = Bh'; -% Row: nvar equations -% Column: 1st lag (with nvar variables) to lags (with nvar variables) + const = k. - -imf = zeros(imstep,nvar*nvar); -% Column: nvar responses to 1st shock, nvar responses to 2nd shock, and so on. -% Row: steps of impulse responses. -M = zeros(nvar*(lags+1),nvar); -% Stack lags M's in the order of, e.g., [Mlags, ..., M2,M1;M0] -M(1:nvar,:) = swish'; -Mtem = M(1:nvar,:); % temporary M. -% first (initial) responses to 1 standard deviation shock. Row: responses; Column: shocks -% * put in the form of "imf" -imf(1,:) = Mtem(:)'; - -t = 1; -ims1 = min([imstep-1 lags]); -while t <= ims1 - Mtem = Ah(:,1:nvar*t)*M(1:nvar*t,:); - % Row: nvar equations, each for the nvar variables at tth lag - M(nvar+1:nvar*(t+1),:)=M(1:nvar*t,:); - M(1:nvar,:) = Mtem; - imf(t+1,:) = Mtem(:)'; - % stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc. - t= t+1; -end - -for t = lags+1:imstep-1 - Mtem = Ah(:,1:nvar*lags)*M(1:nvar*lags,:); - % Row: nvar equations, each for the nvar variables at tth lag - M(nvar+1:nvar*(t+1),:) = M(1:nvar*t,:); - M(1:nvar,:)=Mtem; - imf(t+1,:) = Mtem(:)'; - % stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc. -end +function imf = fn_impulse(Bh,swish,nn) +% Computing impulse functions with +% imf = fn_impulse(Bh,swish,nn) +% imf is in a format that is the SAME as in RATS. +% Column: nvar responses to 1st shock, +% nvar responses to 2nd shock, and so on. +% Row: steps of impulse responses. +%----------------- +% Bh is the estimated reduced form coefficient in the form +% Y(T*nvar) = XB + U, X: T*k (may include all exogenous terms), B: k*nvar. +% The matrix form and dimension are the same as "Bh" from the function "sye.m"; +% Column: 1st lag (with nvar variables) to lags (with nvar variables) + const = k. +% Note: columns correspond to equations. +% swish is the inv(A0) in the structural model y(t)A0 = e(t). +% Note: columns corresponding to equations. +% nn is the numbers of inputs [nvar,lags,# of steps of impulse responses]. +% +% Written by Tao Zha. +% Copyright (c) 1994 by Tao Zha + +% Copyright (C) 1994-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +nvar = nn(1); +lags = nn(2); +imstep = nn(3); % number of steps for impulse responses + +Ah = Bh'; +% Row: nvar equations +% Column: 1st lag (with nvar variables) to lags (with nvar variables) + const = k. + +imf = zeros(imstep,nvar*nvar); +% Column: nvar responses to 1st shock, nvar responses to 2nd shock, and so on. +% Row: steps of impulse responses. +M = zeros(nvar*(lags+1),nvar); +% Stack lags M's in the order of, e.g., [Mlags, ..., M2,M1;M0] +M(1:nvar,:) = swish'; +Mtem = M(1:nvar,:); % temporary M. +% first (initial) responses to 1 standard deviation shock. Row: responses; Column: shocks +% * put in the form of "imf" +imf(1,:) = Mtem(:)'; + +t = 1; +ims1 = min([imstep-1 lags]); +while t <= ims1 + Mtem = Ah(:,1:nvar*t)*M(1:nvar*t,:); + % Row: nvar equations, each for the nvar variables at tth lag + M(nvar+1:nvar*(t+1),:)=M(1:nvar*t,:); + M(1:nvar,:) = Mtem; + imf(t+1,:) = Mtem(:)'; + % stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc. + t= t+1; +end + +for t = lags+1:imstep-1 + Mtem = Ah(:,1:nvar*lags)*M(1:nvar*lags,:); + % Row: nvar equations, each for the nvar variables at tth lag + M(nvar+1:nvar*(t+1),:) = M(1:nvar*t,:); + M(1:nvar,:)=Mtem; + imf(t+1,:) = Mtem(:)'; + % stack imf with each step, Row: 6 var to 1st shock, 6 var to 2nd shock, etc. +end diff --git a/matlab/ms-sbvar/cstz/fn_rlrpostr.m b/matlab/ms-sbvar/cstz/fn_rlrpostr.m index a52653514ca5bb4045880ec726f02224741a8308..f8e5b12a5575402b0ddc336e6b9eb4721b35d6d9 100644 --- a/matlab/ms-sbvar/cstz/fn_rlrpostr.m +++ b/matlab/ms-sbvar/cstz/fn_rlrpostr.m @@ -1,67 +1,67 @@ -function [P,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0invtld,Hpinvtld,Ui,Vi) -% [P,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0tld,Hptld,Ui,Vi) -% -% Exporting random (i.e., random prior) Bayesian posterior matrices with linear restrictions -% See Waggoner and Zha's Gibbs sampling paper -% -% xtx: X'X: k-by-k where k=ncoef -% xty: X'Y: k-by-nvar -% yty: Y'Y: nvar-by-nvar -% Ptld: cell(nvar,1), transformation matrix that affects the (random walk) prior mean of A+ conditional on A0. -% H0invtld: cell(nvar,1), transformed inv covaraince for free parameters in A0(:,i). -% Hpinvtld: cell(nvar,1), transformed inv covaraince for free parameters in A+(:,i); -% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith -% equation contemporaneous restriction matrix where qi is the number of free parameters. -% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector -% of total original parameters and bi is a vector of free parameters. When no -% restrictions are imposed, we have Ui = I. There must be at least one free -% parameter left for the ith equation. Imported from dnrprior.m. -% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith -% equation lagged restriction matrix where k (ncoef) is a total number of RHS variables and -% ri is the number of free parameters. With this transformation, we have fi = Vi*gi -% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a -% vector of free parameters. There must be at least one free parameter left for -% the ith equation. Imported from dnrprior.m. -%----------------- -% P: cell(nvar,1). In each cell, the transformation matrix that affects the posterior mean of A+ conditional on A0. -% In other words, the posterior mean (of g_i) = P{i}*b_i where g_i is a column vector of free parameters -% of A+(:,i)) given b_i (b_i is a column vector of free parameters of A0(:,i)). -% H0inv: cell(nvar,1). Not divided by T yet. In each cell, inverse of posterior covariance matrix H0. -% The exponential term is b_i'*inv(H0)*b_i for the ith equation where b_i = U_i*a0_i. -% It resembles old SpH or Sbd in the exponent term in posterior of A0, but not divided by T yet. -% Hpinv: cell(nvar,1). In each cell, posterior inverse of covariance matrix Hp (A+) for the free parameters -% g_i = V_i*A+(:,i) in the ith equation. -% -% Tao Zha, February 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -nvar = size(yty,1); - -P = cell(nvar,1); % tld: tilda -H0inv = cell(nvar,1); % posterior inv(H0), resemble old SpH, but not divided by T yet. -Hpinv = cell(nvar,1); % posterior inv(Hp). - - -for n=1:nvar % one for each equation - Hpinv{n} = Vi{n}'*xtx*Vi{n} + Hpinvtld{n}; - P1 = Vi{n}'*xty*Ui{n} + Hpinvtld{n}*Ptld{n}; - P{n} = Hpinv{n}\P1; - H0inv{n} = Ui{n}'*yty*Ui{n} + H0invtld{n} + Ptld{n}'*Hpinvtld{n}*Ptld{n} ... - - P1'*P{n}; %P{n} = (Hpinv{n}\P1); -end +function [P,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0invtld,Hpinvtld,Ui,Vi) +% [P,H0inv,Hpinv] = fn_rlrpostr(xtx,xty,yty,Ptld,H0tld,Hptld,Ui,Vi) +% +% Exporting random (i.e., random prior) Bayesian posterior matrices with linear restrictions +% See Waggoner and Zha's Gibbs sampling paper +% +% xtx: X'X: k-by-k where k=ncoef +% xty: X'Y: k-by-nvar +% yty: Y'Y: nvar-by-nvar +% Ptld: cell(nvar,1), transformation matrix that affects the (random walk) prior mean of A+ conditional on A0. +% H0invtld: cell(nvar,1), transformed inv covaraince for free parameters in A0(:,i). +% Hpinvtld: cell(nvar,1), transformed inv covaraince for free parameters in A+(:,i); +% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith +% equation contemporaneous restriction matrix where qi is the number of free parameters. +% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector +% of total original parameters and bi is a vector of free parameters. When no +% restrictions are imposed, we have Ui = I. There must be at least one free +% parameter left for the ith equation. Imported from dnrprior.m. +% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith +% equation lagged restriction matrix where k (ncoef) is a total number of RHS variables and +% ri is the number of free parameters. With this transformation, we have fi = Vi*gi +% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a +% vector of free parameters. There must be at least one free parameter left for +% the ith equation. Imported from dnrprior.m. +%----------------- +% P: cell(nvar,1). In each cell, the transformation matrix that affects the posterior mean of A+ conditional on A0. +% In other words, the posterior mean (of g_i) = P{i}*b_i where g_i is a column vector of free parameters +% of A+(:,i)) given b_i (b_i is a column vector of free parameters of A0(:,i)). +% H0inv: cell(nvar,1). Not divided by T yet. In each cell, inverse of posterior covariance matrix H0. +% The exponential term is b_i'*inv(H0)*b_i for the ith equation where b_i = U_i*a0_i. +% It resembles old SpH or Sbd in the exponent term in posterior of A0, but not divided by T yet. +% Hpinv: cell(nvar,1). In each cell, posterior inverse of covariance matrix Hp (A+) for the free parameters +% g_i = V_i*A+(:,i) in the ith equation. +% +% Tao Zha, February 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +nvar = size(yty,1); + +P = cell(nvar,1); % tld: tilda +H0inv = cell(nvar,1); % posterior inv(H0), resemble old SpH, but not divided by T yet. +Hpinv = cell(nvar,1); % posterior inv(Hp). + + +for n=1:nvar % one for each equation + Hpinv{n} = Vi{n}'*xtx*Vi{n} + Hpinvtld{n}; + P1 = Vi{n}'*xty*Ui{n} + Hpinvtld{n}*Ptld{n}; + P{n} = Hpinv{n}\P1; + H0inv{n} = Ui{n}'*yty*Ui{n} + H0invtld{n} + Ptld{n}'*Hpinvtld{n}*Ptld{n} ... + - P1'*P{n}; %P{n} = (Hpinv{n}\P1); +end diff --git a/matlab/ms-sbvar/cstz/fn_rlrprior.m b/matlab/ms-sbvar/cstz/fn_rlrprior.m index 118b6750962a090c8758d743c08a4c8d5ff78c7f..38600bee9f9a29af7005fadb524c82736d6fa02c 100644 --- a/matlab/ms-sbvar/cstz/fn_rlrprior.m +++ b/matlab/ms-sbvar/cstz/fn_rlrprior.m @@ -1,56 +1,56 @@ -function [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Ui,Vi,Pi,H0multi,Hpmulti,nvar) -% [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Ui,Vi,Pi,H0multi,Hpmulti,nvar) -% -% Exporting random Bayesian prior with linear restrictions -% See Waggoner and Zha's Gibbs sampling paper -% -% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith -% equation contemporaneous restriction matrix where qi is the number of free parameters. -% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector -% of total original parameters and bi is a vector of free parameters. When no -% restrictions are imposed, we have Ui = I. There must be at least one free -% parameter left for the ith equation. Imported from dnrprior.m. -% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith -% equation lagged restriction matrix where k (ncoef) is a total number of RHS variables and -% ri is the number of free parameters. With this transformation, we have fi = Vi*gi -% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a -% vector of free parameters. There must be at least one free parameter left for -% the ith equation. Imported from dnrprior.m. -% Pi: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations -% H0multi: nvar-by-nvar-by-nvar; H0 for different equations under asymmetric prior -% Hpmulti: ncoef-by-ncoef-by-nvar; H+ for different equations under asymmetric prior -% nvar: number of endogenous variables -% -------------------- -% Ptld: cell(nvar,1). The prior mean of g_i is Ptld{i}*b_i; -% H0invtld: cell(nvar,1). Transformed inv covaraince for b_i, the free parameters in A0(:,i); -% Hpinvtld: cell(nvar,1). Transformed inv covaraince for g_i, the free parameters in A+(:,i); -% -% Tao Zha, February 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -Ptld = cell(nvar,1); % tld: tilda -H0invtld = cell(nvar,1); % H0 for different equations under linear restrictions -Hpinvtld = cell(nvar,1); % H+ for different equations under linear restrictions - -for n=1:nvar % one for each equation - Hpinvtld{n} = Vi{n}'*(Hpmulti(:,:,n)\Vi{n}); - Ptld{n} = (Hpinvtld{n}\Vi{n}')*(Hpmulti(:,:,n)\Pi)*Ui{n}; - H0invtld{n} = Ui{n}'*(H0multi(:,:,n)\Ui{n}) + Ui{n}'*Pi'*(Hpmulti(:,:,n)\Pi)*Ui{n} ... - - Ptld{n}'*Hpinvtld{n}*Ptld{n}; -end +function [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Ui,Vi,Pi,H0multi,Hpmulti,nvar) +% [Ptld,H0invtld,Hpinvtld] = fn_rlrprior(Ui,Vi,Pi,H0multi,Hpmulti,nvar) +% +% Exporting random Bayesian prior with linear restrictions +% See Waggoner and Zha's Gibbs sampling paper +% +% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith +% equation contemporaneous restriction matrix where qi is the number of free parameters. +% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector +% of total original parameters and bi is a vector of free parameters. When no +% restrictions are imposed, we have Ui = I. There must be at least one free +% parameter left for the ith equation. Imported from dnrprior.m. +% Vi: nvar-by-1 cell. In each cell, k-by-ri orthonormal basis for the null of the ith +% equation lagged restriction matrix where k (ncoef) is a total number of RHS variables and +% ri is the number of free parameters. With this transformation, we have fi = Vi*gi +% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a +% vector of free parameters. There must be at least one free parameter left for +% the ith equation. Imported from dnrprior.m. +% Pi: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations +% H0multi: nvar-by-nvar-by-nvar; H0 for different equations under asymmetric prior +% Hpmulti: ncoef-by-ncoef-by-nvar; H+ for different equations under asymmetric prior +% nvar: number of endogenous variables +% -------------------- +% Ptld: cell(nvar,1). The prior mean of g_i is Ptld{i}*b_i; +% H0invtld: cell(nvar,1). Transformed inv covaraince for b_i, the free parameters in A0(:,i); +% Hpinvtld: cell(nvar,1). Transformed inv covaraince for g_i, the free parameters in A+(:,i); +% +% Tao Zha, February 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +Ptld = cell(nvar,1); % tld: tilda +H0invtld = cell(nvar,1); % H0 for different equations under linear restrictions +Hpinvtld = cell(nvar,1); % H+ for different equations under linear restrictions + +for n=1:nvar % one for each equation + Hpinvtld{n} = Vi{n}'*(Hpmulti(:,:,n)\Vi{n}); + Ptld{n} = (Hpinvtld{n}\Vi{n}')*(Hpmulti(:,:,n)\Pi)*Ui{n}; + H0invtld{n} = Ui{n}'*(H0multi(:,:,n)\Ui{n}) + Ui{n}'*Pi'*(Hpmulti(:,:,n)\Pi)*Ui{n} ... + - Ptld{n}'*Hpinvtld{n}*Ptld{n}; +end diff --git a/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m b/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m index 7986d48376fdd03287c3dc38709582e27ac3a899..394d27b90ffe3faa779a82fd9eba06c86b7fa8f0 100644 --- a/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m +++ b/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs.m @@ -1,295 +1,295 @@ -function [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] ... - = fn_rnrprior_covres_dobs(nvar,q_m,lags,xdgel,mu,indxDummy,hpmsmd,indxmsmdeqn,nexo,asym0,asymp) -% Differs from fn_rnrprior_covres_dobs_tv(): no linear restrictions (Ui and Vi) have applied yet to this function, but -% linear restrictions are incorported in fn_rnrprior_covres_dobs_tv(). -% -% Only works for the nexo=1 (constant term) case. To extend this to other exogenous variables, see fn_dataxy.m. 01/14/03. -% Differs from fn_rnrprior_covres.m in that dummy observations are included as part of the explicit prior. See Forcast II, pp.68-69b. -% More general than fn_rnrprior.m because when hpmsmd=0, fn_rnrprior_covres() is the same as fn_rnrprior(). -% Allows for prior covariances for the MS and MD equations to achieve liquidity effects. -% Exports random Bayesian prior of Sims and Zha with asymmetric rior (but no linear restrictions yet) -% See Waggoner and Zha's Gibbs sampling paper and TVBVAR NOTES p. 71k.0. -% -% nvar: number of endogenous variables -% q_m: quarter or month -% lags: the maximum length of lag -% xdgel: T*nvar endogenous-variable matrix of raw or original data (no manipulation involved) with sample size including lags. -% Order of columns: (1) nvar endogenous variables; (2) constants will be automatically put in the last column. -% Used only to get variances of residuals for mu(1)-mu(5) and for dummy observations mu(5) and mu(6). -% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where -% mu(5) and mu(6) are NOT used here. See fn_dataxy.m for using mu(5) and mu(6). -% mu(1): overall tightness and also for A0; (0.57) -% mu(2): relative tightness for A+; (0.13) -% mu(3): relative tightness for the constant term; (0.1). NOTE: for other -% exogenous terms, the variance of each exogenous term must be taken into -% acount to eliminate the scaling factor. -% mu(4): tightness on lag decay; (1) -% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5) -% mu(6): weight on single dummy initial observation including constant -% (cointegration, unit roots, and stationarity); (5) -% indxDummy: 1: uses dummy observations to form part of an explicit prior; 0: no dummy observations as part of the prior. -% hpmsmd: 2-by-1 hyperparameters with -1<h1=hpmsmd(1)<=0 for the MS equation and 0<=h2=hpmsmd(2)<1 the MD equation. Consider a1*R + a2*M. -% The term h1*var(a1)*var(a2) is the prior covariance of a1 and a2 for MS, equivalent to penalizing the same sign of a1 and a2. -% The term h2*var(a1)*var(a2) is the prior covariance of a1 and a2 for MD, equivalent to penalizing opposite signs of a1 and a2. -% This will give us a liquidity effect. -% indxmsmdeqn: 4-by-1 index for the locations of the MS and MD equation and for the locations of M and R. -% indxmsmdeqn(1) for MS and indxmsmdeqn(2) for MD. -% indxmsmdeqn(3) for M and indxmsmdeqn(4) for R. -% nexo: number of exogenous variables (if not specified, nexo=1 (constant) by default). -% The constant term is always put to the last of all endogenous and exogenous variables. -% asym0: nvar-by-nvar asymmetric prior on A0. Column -- equation. -% If ones(nvar,nvar), symmetric prior; if not, relative (asymmetric) tightness on A0. -% asymp: ncoef-1-by-nvar asymmetric prior on A+ bar constant. Column -- equation. -% If ones(ncoef-1,nvar), symmetric prior; if not, relative (asymmetric) tightness on A+. -% -------------------- -% Pi: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations -% H0multi: nvar-by-nvar-by-nvar; H0 for different equations under asymmetric prior -% Hpmulti: ncoef-by-ncoef-by-nvar; H+ for different equations under asymmetric prior -% H0invmulti: nvar-by-nvar-by-nvar; inv(H0) for different equations under asymmetric prior -% Hpinvmulti: ncoef-by-ncoef-by-nvar; inv(H+) for different equations under asymmetric prior -% -% Tao Zha, February 2000. Revised, September 2000, February, May 2003. - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if (nargin<=8), nexo=1; end -ncoef = nvar*lags+nexo; % number of coefficients in *each* equation, RHS coefficients only. - -H0multi=zeros(nvar,nvar,nvar); % H0 for different equations under asymmetric prior -Hpmulti=zeros(ncoef,ncoef,nvar); % H+ for different equations under asymmetric prior -H0invmulti=zeros(nvar,nvar,nvar); % inv(H0) for different equations under asymmetric prior -Hpinvmulti=zeros(ncoef,ncoef,nvar); % inv(H+) for different equations under asymmetric prior - -%*** Constructing Pi for the ith equation under the random walk assumption -Pi = zeros(ncoef,nvar); % same for all equations -Pi(1:nvar,1:nvar) = eye(nvar); % random walk - -% -%@@@ Prepared for Bayesian prior -% -% -% ** monthly lag decay in order to match quarterly decay: a*exp(bl) where -% ** l is the monthly lag. Suppose quarterly decay is 1/x where x=1,2,3,4. -% ** Let the decay of l1 (a*exp(b*l1)) match that of x1 (say, beginning: 1/1) -% ** and the decay of l2 (a*exp(b*l2)) match that of x2 (say, end: 1/5), -% ** we can solve for a and b which are -% ** b = (log_x1-log_x2)/(l1-l2), and a = x1*exp(-b*l1). -if q_m==12 - l1 = 1; % 1st month == 1st quarter - xx1 = 1; % 1st quarter - l2 = lags; % last month - xx2 = 1/((ceil(lags/3))^mu(4)); % last quarter - %xx2 = 1/6; % last quarter - % 3rd quarter: i.e., we intend to let decay of the 6th month match - % that of the 3rd quarter, so that the 6th month decays a little - % faster than the second quarter which is 1/2. - if lags==1 - b = 0; - else - b = (log(xx1)-log(xx2))/(l1-l2); - end - a = xx1*exp(-b*l1); -end - - - -% -% *** specify the prior for each equation separately, SZ method, -% ** get the residuals from univariate regressions. -% -sgh = zeros(nvar,1); % square root -sgsh = sgh; % square -nSample=size(xdgel,1); % sample size-lags -yu = xdgel; -C = ones(nSample,1); -for k=1:nvar - [Bk,ek,junk1,junk2,junk3,junk4] = sye([yu(:,k) C],lags); - clear Bk junk1 junk2 junk3 junk4; - sgsh(k) = ek'*ek/(nSample-lags); - sgh(k) = sqrt(sgsh(k)); -end -% ** prior variance for A0(:,1), same for all equations!!! -sg0bid = zeros(nvar,1); % Sigma0_bar diagonal only for the ith equation -for j=1:nvar - sg0bid(j) = 1/sgsh(j); % sgsh = sigmai^2 -end -% ** prior variance for lagged and exogeous variables, same for all equations -sgpbid = zeros(ncoef,1); % Sigma_plus_bar, diagonal, for the ith equation -for i = 1:lags - if (q_m==12) - lagdecay = a*exp(b*i*mu(4)); - end - % - for j = 1:nvar - if (q_m==12) - % exponential decay to match quarterly decay - sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation - elseif (q_m==4) - sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation - 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 - - +function [Pi,H0multi,Hpmulti,H0invmulti,Hpinvmulti] ... + = fn_rnrprior_covres_dobs(nvar,q_m,lags,xdgel,mu,indxDummy,hpmsmd,indxmsmdeqn,nexo,asym0,asymp) +% Differs from fn_rnrprior_covres_dobs_tv(): no linear restrictions (Ui and Vi) have applied yet to this function, but +% linear restrictions are incorported in fn_rnrprior_covres_dobs_tv(). +% +% Only works for the nexo=1 (constant term) case. To extend this to other exogenous variables, see fn_dataxy.m. 01/14/03. +% Differs from fn_rnrprior_covres.m in that dummy observations are included as part of the explicit prior. See Forcast II, pp.68-69b. +% More general than fn_rnrprior.m because when hpmsmd=0, fn_rnrprior_covres() is the same as fn_rnrprior(). +% Allows for prior covariances for the MS and MD equations to achieve liquidity effects. +% Exports random Bayesian prior of Sims and Zha with asymmetric rior (but no linear restrictions yet) +% See Waggoner and Zha's Gibbs sampling paper and TVBVAR NOTES p. 71k.0. +% +% nvar: number of endogenous variables +% q_m: quarter or month +% lags: the maximum length of lag +% xdgel: T*nvar endogenous-variable matrix of raw or original data (no manipulation involved) with sample size including lags. +% Order of columns: (1) nvar endogenous variables; (2) constants will be automatically put in the last column. +% Used only to get variances of residuals for mu(1)-mu(5) and for dummy observations mu(5) and mu(6). +% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where +% mu(5) and mu(6) are NOT used here. See fn_dataxy.m for using mu(5) and mu(6). +% mu(1): overall tightness and also for A0; (0.57) +% mu(2): relative tightness for A+; (0.13) +% mu(3): relative tightness for the constant term; (0.1). NOTE: for other +% exogenous terms, the variance of each exogenous term must be taken into +% acount to eliminate the scaling factor. +% mu(4): tightness on lag decay; (1) +% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5) +% mu(6): weight on single dummy initial observation including constant +% (cointegration, unit roots, and stationarity); (5) +% indxDummy: 1: uses dummy observations to form part of an explicit prior; 0: no dummy observations as part of the prior. +% hpmsmd: 2-by-1 hyperparameters with -1<h1=hpmsmd(1)<=0 for the MS equation and 0<=h2=hpmsmd(2)<1 the MD equation. Consider a1*R + a2*M. +% The term h1*var(a1)*var(a2) is the prior covariance of a1 and a2 for MS, equivalent to penalizing the same sign of a1 and a2. +% The term h2*var(a1)*var(a2) is the prior covariance of a1 and a2 for MD, equivalent to penalizing opposite signs of a1 and a2. +% This will give us a liquidity effect. +% indxmsmdeqn: 4-by-1 index for the locations of the MS and MD equation and for the locations of M and R. +% indxmsmdeqn(1) for MS and indxmsmdeqn(2) for MD. +% indxmsmdeqn(3) for M and indxmsmdeqn(4) for R. +% nexo: number of exogenous variables (if not specified, nexo=1 (constant) by default). +% The constant term is always put to the last of all endogenous and exogenous variables. +% asym0: nvar-by-nvar asymmetric prior on A0. Column -- equation. +% If ones(nvar,nvar), symmetric prior; if not, relative (asymmetric) tightness on A0. +% asymp: ncoef-1-by-nvar asymmetric prior on A+ bar constant. Column -- equation. +% If ones(ncoef-1,nvar), symmetric prior; if not, relative (asymmetric) tightness on A+. +% -------------------- +% Pi: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations +% H0multi: nvar-by-nvar-by-nvar; H0 for different equations under asymmetric prior +% Hpmulti: ncoef-by-ncoef-by-nvar; H+ for different equations under asymmetric prior +% H0invmulti: nvar-by-nvar-by-nvar; inv(H0) for different equations under asymmetric prior +% Hpinvmulti: ncoef-by-ncoef-by-nvar; inv(H+) for different equations under asymmetric prior +% +% Tao Zha, February 2000. Revised, September 2000, February, May 2003. + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if (nargin<=8), nexo=1; end +ncoef = nvar*lags+nexo; % number of coefficients in *each* equation, RHS coefficients only. + +H0multi=zeros(nvar,nvar,nvar); % H0 for different equations under asymmetric prior +Hpmulti=zeros(ncoef,ncoef,nvar); % H+ for different equations under asymmetric prior +H0invmulti=zeros(nvar,nvar,nvar); % inv(H0) for different equations under asymmetric prior +Hpinvmulti=zeros(ncoef,ncoef,nvar); % inv(H+) for different equations under asymmetric prior + +%*** Constructing Pi for the ith equation under the random walk assumption +Pi = zeros(ncoef,nvar); % same for all equations +Pi(1:nvar,1:nvar) = eye(nvar); % random walk + +% +%@@@ Prepared for Bayesian prior +% +% +% ** monthly lag decay in order to match quarterly decay: a*exp(bl) where +% ** l is the monthly lag. Suppose quarterly decay is 1/x where x=1,2,3,4. +% ** Let the decay of l1 (a*exp(b*l1)) match that of x1 (say, beginning: 1/1) +% ** and the decay of l2 (a*exp(b*l2)) match that of x2 (say, end: 1/5), +% ** we can solve for a and b which are +% ** b = (log_x1-log_x2)/(l1-l2), and a = x1*exp(-b*l1). +if q_m==12 + l1 = 1; % 1st month == 1st quarter + xx1 = 1; % 1st quarter + l2 = lags; % last month + xx2 = 1/((ceil(lags/3))^mu(4)); % last quarter + %xx2 = 1/6; % last quarter + % 3rd quarter: i.e., we intend to let decay of the 6th month match + % that of the 3rd quarter, so that the 6th month decays a little + % faster than the second quarter which is 1/2. + if lags==1 + b = 0; + else + b = (log(xx1)-log(xx2))/(l1-l2); + end + a = xx1*exp(-b*l1); +end + + + +% +% *** specify the prior for each equation separately, SZ method, +% ** get the residuals from univariate regressions. +% +sgh = zeros(nvar,1); % square root +sgsh = sgh; % square +nSample=size(xdgel,1); % sample size-lags +yu = xdgel; +C = ones(nSample,1); +for k=1:nvar + [Bk,ek,junk1,junk2,junk3,junk4] = sye([yu(:,k) C],lags); + clear Bk junk1 junk2 junk3 junk4; + sgsh(k) = ek'*ek/(nSample-lags); + sgh(k) = sqrt(sgsh(k)); +end +% ** prior variance for A0(:,1), same for all equations!!! +sg0bid = zeros(nvar,1); % Sigma0_bar diagonal only for the ith equation +for j=1:nvar + sg0bid(j) = 1/sgsh(j); % sgsh = sigmai^2 +end +% ** prior variance for lagged and exogeous variables, same for all equations +sgpbid = zeros(ncoef,1); % Sigma_plus_bar, diagonal, for the ith equation +for i = 1:lags + if (q_m==12) + lagdecay = a*exp(b*i*mu(4)); + end + % + for j = 1:nvar + if (q_m==12) + % exponential decay to match quarterly decay + sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation + elseif (q_m==4) + sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation + 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 index 66c99a93562c3f42fe2b20c7a3e7e2547b467602..e3d10d4094baedd3fb25aafb0e43d6b43fa4d718 100644 --- a/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs_tv2.m +++ b/matlab/ms-sbvar/cstz/fn_rnrprior_covres_dobs_tv2.m @@ -1,325 +1,325 @@ -function [Pi_bar,H0tldcell_inv,Hptldcell_inv] ... - = fn_rnrprior_covres_dobs_tv2(nvar,nStates,indxScaleStates,q_m,lags,xdgel,mu,indxDummy,Ui,Vi,hpmsmd,indxmsmdeqn,nexo,asym0,asymp) -% Differs from fn_rnrprior_covres_dobs(): linear restrictions (Ui and Vi) have been incorported in fn_rnrprior_covres_dobs_tv?(). -% Differs from fn_rnrprior_covres_dobs_tv(): allows an option to scale up the prior variance by nStates or not scale at all, -% so that the prior value is the same as the constant VAR when the parameters in all states are the same. -% -% Only works for the nexo=1 (constant term) case. To extend this to other exogenous variables, see fn_dataxy.m. 01/14/03. -% Differs from fn_rnrprior_covres_tv.m in that dummy observations are included as part of the explicit prior. See Forcast II, pp.68-69b. -% Exports random Bayesian prior of Sims and Zha with asymmetric rior with linear restrictions already applied -% and with dummy observations (i.e., mu(5) and mu(6)) used as part of an explicit prior. -% This function allows for prior covariances for the MS and MD equations to achieve liquidity effects. -% See Waggoner and Zha's Gibbs sampling paper and TVBVAR NOTES pp. 71k.0 and 50-61. -% -% nvar: number of endogenous variables -% nStates: Number of states. -% indxScaleStates: if 0, no scale adjustment in the prior variance for the number of states in the function fn_rnrprior_covres_dobs_tv2(); -% if 1: allows a scale adjustment, marking the prior variance bigger by the number of states. -% q_m: quarter or month -% lags: the maximum length of lag -% xdgel: T*nvar endogenous-variable matrix of raw or original data (no manipulation involved) with sample size including lags. -% Order of columns: (1) nvar endogenous variables; (2) constants will be automatically put in the last column. -% Used only to get variances of residuals for mu(1)-mu(5) and for dummy observations mu(5) and mu(6). -% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where -% mu(5) and mu(6) are NOT used here. See fn_dataxy.m for using mu(5) and mu(6). -% mu(1): overall tightness and also for A0; (0.57) -% mu(2): relative tightness for A+; (0.13) -% mu(3): relative tightness for the constant term; (0.1). NOTE: for other -% exogenous terms, the variance of each exogenous term must be taken into -% acount to eliminate the scaling factor. -% mu(4): tightness on lag decay; (1) -% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5) -% mu(6): weight on single dummy initial observation including constant -% (cointegration, unit roots, and stationarity); (5) -% NOTE: for this function, mu(5) and mu(6) are not used. See fn_dataxy.m for using mu(5) and mu(6). -% indxDummy: 1: uses dummy observations to form part of an explicit prior; 0: no dummy observations as part of the prior. -% Ui: nvar-by-1 cell. In each cell, nvar-by-qi*si orthonormal basis for the null of the ith -% equation contemporaneous restriction matrix where qi is the number of free parameters -% within the state and si is the number of free states. -% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector -% of total original parameters and bi is a vector of free parameters. When no -% restrictions are imposed, we have Ui = I. There must be at least one free -% parameter left for the ith equation in the order of [a_i for 1st state, ..., a_i for last state]. -% Vi: nvar-by-1 cell. In each cell, k-by-ri*ti orthonormal basis for the null of the ith -% equation lagged restriction matrix where k is a total of exogenous variables and -% ri is the number of free parameters within the state and ti is the number of free states. -% With this transformation, we have fi = Vi*gi -% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a -% vector of free parameters. The ith equation is in the order of [nvar variables -% for 1st lag and 1st state, ..., nvar variables for last lag and 1st state, const for 1st state, nvar -% variables for 1st lag and 2nd state, nvar variables for last lag and 2nd state, const for 2nd state, and so on]. -% hpmsmd: 2-by-1 hyperparameters with -1<h1=hpmsmd(1)<=0 for the MS equation and 0<=h2=hpmsmd(2)<1 the MD equation. Consider a1*R + a2*M. -% The term h1*var(a1)*var(a2) is the prior covariance of a1 and a2 for MS, equivalent to penalizing the same sign of a1 and a2. -% The term h2*var(a1)*var(a2) is the prior covariance of a1 and a2 for MD, equivalent to penalizing opposite signs of a1 and a2. -% This will give us a liquidity effect. If hpmsmd=0, no such restrictions will be imposed. -% indxmsmdeqn: 4-by-1 index for the locations of the MS and MD equation and for the locations of M and R. -% indxmsmdeqn(1) for MS and indxmsmdeqn(2) for MD. -% indxmsmdeqn(3) for M and indxmsmdeqn(4) for R. -% nexo: number of exogenous variables (if not specified, nexo=1 (constant) by default). -% The constant term is always put to the last of all endogenous and exogenous variables. -% asym0: nvar-by-nvar asymmetric prior on A0. Column -- equation. -% If ones(nvar,nvar), symmetric prior; if not, relative (asymmetric) tightness on A0. -% asymp: ncoef-1-by-nvar asymmetric prior on A+ bar constant. Column -- equation. -% If ones(ncoef-1,nvar), symmetric prior; if not, relative (asymmetric) tightness on A+. -% -------------------- -% Pi_bar: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations -% H0tldcell_inv: cell(nvar,1). The ith cell represents the ith equation, where the dim is -% qi*si-by-qi*si. The inverse of H0tld on p.60. -% Hptldcell_inv: cell(nvar,1). The ith cell represents the ith equation, where the dim is -% ri*ti-by-ri*ti.The inverse of Hptld on p.60. -% -% Differs from fn_rnrprior_covres_dobs(): linear restrictions (Ui and Vi) have been incorported in fn_rnrprior_covres_dobs_tv?(). -% Differs from fn_rnrprior_covres_dobs_tv(): allows an option to scale up the prior variance by nStates or not scale at all. -% so that the prior value is the same as the constant VAR when the parameters in all states are the same. -% Tao Zha, February 2000. Revised, September 2000, 2001, February, May 2003, May 2004. - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if nargin<=12, nexo=1; end % <<>>1 -ncoef = nvar*lags+nexo; % Number of coefficients in *each* equation for each state, RHS coefficients only. -ncoefsts = nStates*ncoef; % Number of coefficients in *each* equation in all states, RHS coefficients only. - -H0tldcell_inv=cell(nvar,1); % inv(H0tilde) for different equations under asymmetric prior. -Hptldcell_inv=cell(nvar,1); % inv(H+tilde) for different equations under asymmetric prior. - -%*** Constructing Pi_bar for the ith equation under the random walk assumption -Pi_bar = zeros(ncoef,nvar); % same for all equations -Pi_bar(1:nvar,1:nvar) = eye(nvar); % random walk - -% -%@@@ Prepared for Bayesian prior -% -% -% ** monthly lag decay in order to match quarterly decay: a*exp(bl) where -% ** l is the monthly lag. Suppose quarterly decay is 1/x where x=1,2,3,4. -% ** Let the decay of l1 (a*exp(b*l1)) match that of x1 (say, beginning: 1/1) -% ** and the decay of l2 (a*exp(b*l2)) match that of x2 (say, end: 1/5), -% ** we can solve for a and b which are -% ** b = (log_x1-log_x2)/(l1-l2), and a = x1*exp(-b*l1). -if q_m==12 - l1 = 1; % 1st month == 1st quarter - xx1 = 1; % 1st quarter - l2 = lags; % last month - xx2 = 1/((ceil(lags/3))^mu(4)); % last quarter - %xx2 = 1/6; % last quarter - % 3rd quarter: i.e., we intend to let decay of the 6th month match - % that of the 3rd quarter, so that the 6th month decays a little - % faster than the second quarter which is 1/2. - if lags==1 - b = 0; - else - b = (log(xx1)-log(xx2))/(l1-l2); - end - a = xx1*exp(-b*l1); -end - - - -% -% *** specify the prior for each equation separately, SZ method, -% ** get the residuals from univariate regressions. -% -sgh = zeros(nvar,1); % square root -sgsh = sgh; % square -nSample=size(xdgel,1); % sample size-lags -yu = xdgel; -C = ones(nSample,1); -for k=1:nvar - [Bk,ek,junk1,junk2,junk3,junk4] = sye([yu(:,k) C],lags); - clear Bk junk1 junk2 junk3 junk4; - sgsh(k) = ek'*ek/(nSample-lags); - sgh(k) = sqrt(sgsh(k)); -end -% ** prior variance for A0(:,1), same for all equations!!! -sg0bid = zeros(nvar,1); % Sigma0_bar diagonal only for the ith equation -for j=1:nvar - sg0bid(j) = 1/sgsh(j); % sgsh = sigmai^2 -end -% ** prior variance for lagged and exogeous variables, same for all equations -sgpbid = zeros(ncoef,1); % Sigma_plus_bar, diagonal, for the ith equation -for i = 1:lags - if (q_m==12) - lagdecay = a*exp(b*i*mu(4)); - end - % - for j = 1:nvar - if (q_m==12) - % exponential decay to match quarterly decay - sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation - elseif (q_m==4) - sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation - 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 - - +function [Pi_bar,H0tldcell_inv,Hptldcell_inv] ... + = fn_rnrprior_covres_dobs_tv2(nvar,nStates,indxScaleStates,q_m,lags,xdgel,mu,indxDummy,Ui,Vi,hpmsmd,indxmsmdeqn,nexo,asym0,asymp) +% Differs from fn_rnrprior_covres_dobs(): linear restrictions (Ui and Vi) have been incorported in fn_rnrprior_covres_dobs_tv?(). +% Differs from fn_rnrprior_covres_dobs_tv(): allows an option to scale up the prior variance by nStates or not scale at all, +% so that the prior value is the same as the constant VAR when the parameters in all states are the same. +% +% Only works for the nexo=1 (constant term) case. To extend this to other exogenous variables, see fn_dataxy.m. 01/14/03. +% Differs from fn_rnrprior_covres_tv.m in that dummy observations are included as part of the explicit prior. See Forcast II, pp.68-69b. +% Exports random Bayesian prior of Sims and Zha with asymmetric rior with linear restrictions already applied +% and with dummy observations (i.e., mu(5) and mu(6)) used as part of an explicit prior. +% This function allows for prior covariances for the MS and MD equations to achieve liquidity effects. +% See Waggoner and Zha's Gibbs sampling paper and TVBVAR NOTES pp. 71k.0 and 50-61. +% +% nvar: number of endogenous variables +% nStates: Number of states. +% indxScaleStates: if 0, no scale adjustment in the prior variance for the number of states in the function fn_rnrprior_covres_dobs_tv2(); +% if 1: allows a scale adjustment, marking the prior variance bigger by the number of states. +% q_m: quarter or month +% lags: the maximum length of lag +% xdgel: T*nvar endogenous-variable matrix of raw or original data (no manipulation involved) with sample size including lags. +% Order of columns: (1) nvar endogenous variables; (2) constants will be automatically put in the last column. +% Used only to get variances of residuals for mu(1)-mu(5) and for dummy observations mu(5) and mu(6). +% mu: 6-by-1 vector of hyperparameters (the following numbers for Atlanta Fed's forecast), where +% mu(5) and mu(6) are NOT used here. See fn_dataxy.m for using mu(5) and mu(6). +% mu(1): overall tightness and also for A0; (0.57) +% mu(2): relative tightness for A+; (0.13) +% mu(3): relative tightness for the constant term; (0.1). NOTE: for other +% exogenous terms, the variance of each exogenous term must be taken into +% acount to eliminate the scaling factor. +% mu(4): tightness on lag decay; (1) +% mu(5): weight on nvar sums of coeffs dummy observations (unit roots); (5) +% mu(6): weight on single dummy initial observation including constant +% (cointegration, unit roots, and stationarity); (5) +% NOTE: for this function, mu(5) and mu(6) are not used. See fn_dataxy.m for using mu(5) and mu(6). +% indxDummy: 1: uses dummy observations to form part of an explicit prior; 0: no dummy observations as part of the prior. +% Ui: nvar-by-1 cell. In each cell, nvar-by-qi*si orthonormal basis for the null of the ith +% equation contemporaneous restriction matrix where qi is the number of free parameters +% within the state and si is the number of free states. +% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector +% of total original parameters and bi is a vector of free parameters. When no +% restrictions are imposed, we have Ui = I. There must be at least one free +% parameter left for the ith equation in the order of [a_i for 1st state, ..., a_i for last state]. +% Vi: nvar-by-1 cell. In each cell, k-by-ri*ti orthonormal basis for the null of the ith +% equation lagged restriction matrix where k is a total of exogenous variables and +% ri is the number of free parameters within the state and ti is the number of free states. +% With this transformation, we have fi = Vi*gi +% or Vi'*fi = gi where fi is a vector of total original parameters and gi is a +% vector of free parameters. The ith equation is in the order of [nvar variables +% for 1st lag and 1st state, ..., nvar variables for last lag and 1st state, const for 1st state, nvar +% variables for 1st lag and 2nd state, nvar variables for last lag and 2nd state, const for 2nd state, and so on]. +% hpmsmd: 2-by-1 hyperparameters with -1<h1=hpmsmd(1)<=0 for the MS equation and 0<=h2=hpmsmd(2)<1 the MD equation. Consider a1*R + a2*M. +% The term h1*var(a1)*var(a2) is the prior covariance of a1 and a2 for MS, equivalent to penalizing the same sign of a1 and a2. +% The term h2*var(a1)*var(a2) is the prior covariance of a1 and a2 for MD, equivalent to penalizing opposite signs of a1 and a2. +% This will give us a liquidity effect. If hpmsmd=0, no such restrictions will be imposed. +% indxmsmdeqn: 4-by-1 index for the locations of the MS and MD equation and for the locations of M and R. +% indxmsmdeqn(1) for MS and indxmsmdeqn(2) for MD. +% indxmsmdeqn(3) for M and indxmsmdeqn(4) for R. +% nexo: number of exogenous variables (if not specified, nexo=1 (constant) by default). +% The constant term is always put to the last of all endogenous and exogenous variables. +% asym0: nvar-by-nvar asymmetric prior on A0. Column -- equation. +% If ones(nvar,nvar), symmetric prior; if not, relative (asymmetric) tightness on A0. +% asymp: ncoef-1-by-nvar asymmetric prior on A+ bar constant. Column -- equation. +% If ones(ncoef-1,nvar), symmetric prior; if not, relative (asymmetric) tightness on A+. +% -------------------- +% Pi_bar: ncoef-by-nvar matrix for the ith equation under random walk. Same for all equations +% H0tldcell_inv: cell(nvar,1). The ith cell represents the ith equation, where the dim is +% qi*si-by-qi*si. The inverse of H0tld on p.60. +% Hptldcell_inv: cell(nvar,1). The ith cell represents the ith equation, where the dim is +% ri*ti-by-ri*ti.The inverse of Hptld on p.60. +% +% Differs from fn_rnrprior_covres_dobs(): linear restrictions (Ui and Vi) have been incorported in fn_rnrprior_covres_dobs_tv?(). +% Differs from fn_rnrprior_covres_dobs_tv(): allows an option to scale up the prior variance by nStates or not scale at all. +% so that the prior value is the same as the constant VAR when the parameters in all states are the same. +% Tao Zha, February 2000. Revised, September 2000, 2001, February, May 2003, May 2004. + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if nargin<=12, nexo=1; end % <<>>1 +ncoef = nvar*lags+nexo; % Number of coefficients in *each* equation for each state, RHS coefficients only. +ncoefsts = nStates*ncoef; % Number of coefficients in *each* equation in all states, RHS coefficients only. + +H0tldcell_inv=cell(nvar,1); % inv(H0tilde) for different equations under asymmetric prior. +Hptldcell_inv=cell(nvar,1); % inv(H+tilde) for different equations under asymmetric prior. + +%*** Constructing Pi_bar for the ith equation under the random walk assumption +Pi_bar = zeros(ncoef,nvar); % same for all equations +Pi_bar(1:nvar,1:nvar) = eye(nvar); % random walk + +% +%@@@ Prepared for Bayesian prior +% +% +% ** monthly lag decay in order to match quarterly decay: a*exp(bl) where +% ** l is the monthly lag. Suppose quarterly decay is 1/x where x=1,2,3,4. +% ** Let the decay of l1 (a*exp(b*l1)) match that of x1 (say, beginning: 1/1) +% ** and the decay of l2 (a*exp(b*l2)) match that of x2 (say, end: 1/5), +% ** we can solve for a and b which are +% ** b = (log_x1-log_x2)/(l1-l2), and a = x1*exp(-b*l1). +if q_m==12 + l1 = 1; % 1st month == 1st quarter + xx1 = 1; % 1st quarter + l2 = lags; % last month + xx2 = 1/((ceil(lags/3))^mu(4)); % last quarter + %xx2 = 1/6; % last quarter + % 3rd quarter: i.e., we intend to let decay of the 6th month match + % that of the 3rd quarter, so that the 6th month decays a little + % faster than the second quarter which is 1/2. + if lags==1 + b = 0; + else + b = (log(xx1)-log(xx2))/(l1-l2); + end + a = xx1*exp(-b*l1); +end + + + +% +% *** specify the prior for each equation separately, SZ method, +% ** get the residuals from univariate regressions. +% +sgh = zeros(nvar,1); % square root +sgsh = sgh; % square +nSample=size(xdgel,1); % sample size-lags +yu = xdgel; +C = ones(nSample,1); +for k=1:nvar + [Bk,ek,junk1,junk2,junk3,junk4] = sye([yu(:,k) C],lags); + clear Bk junk1 junk2 junk3 junk4; + sgsh(k) = ek'*ek/(nSample-lags); + sgh(k) = sqrt(sgsh(k)); +end +% ** prior variance for A0(:,1), same for all equations!!! +sg0bid = zeros(nvar,1); % Sigma0_bar diagonal only for the ith equation +for j=1:nvar + sg0bid(j) = 1/sgsh(j); % sgsh = sigmai^2 +end +% ** prior variance for lagged and exogeous variables, same for all equations +sgpbid = zeros(ncoef,1); % Sigma_plus_bar, diagonal, for the ith equation +for i = 1:lags + if (q_m==12) + lagdecay = a*exp(b*i*mu(4)); + end + % + for j = 1:nvar + if (q_m==12) + % exponential decay to match quarterly decay + sgpbid((i-1)*nvar+j) = lagdecay^2/sgsh(j); % ith equation + elseif (q_m==4) + sgpbid((i-1)*nvar+j) = (1/i^mu(4))^2/sgsh(j); % ith equation + 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 index 34da07909733295beeefffb858b494e244c28648..c10e2eee513a5d7163c4b19566f7d0b92776425c 100644 --- a/matlab/ms-sbvar/cstz/fn_tran_a2b.m +++ b/matlab/ms-sbvar/cstz/fn_tran_a2b.m @@ -1,42 +1,42 @@ -function b = fn_tran_a2b(A0,Ui,nvar,n0) -% b = fn_tran_a2b(A0,Ui,nvar,n0) -% Transform A0 to free parameters b's. Note: columns correspond to equations -% See Waggoner and Zha's ``A Gibbs sampler for structural VARs'' -% -% A0: nvar-by-nvar, contempareous matrix (columns correspond to equations) -% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith -% equation contemporaneous restriction matrix where qi is the number of free parameters. -% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector -% of total original parameters and bi is a vector of free parameters. When no -% restrictions are imposed, we have Ui = I. There must be at least one free -% parameter left for the ith equation. -% nvar: number of endogeous variables -% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation -%---------------- -% b: sum(n0)-by-1 vector of A0 free parameters -% -% Tao Zha, February 2000. Revised, August 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -n0=n0(:); -n0cum = [0; cumsum(n0)]; -b=zeros(n0cum(end),1); -for kj = 1:nvar - b(n0cum(kj)+1:n0cum(kj+1))=Ui{kj}'*A0(:,kj); -end +function b = fn_tran_a2b(A0,Ui,nvar,n0) +% b = fn_tran_a2b(A0,Ui,nvar,n0) +% Transform A0 to free parameters b's. Note: columns correspond to equations +% See Waggoner and Zha's ``A Gibbs sampler for structural VARs'' +% +% A0: nvar-by-nvar, contempareous matrix (columns correspond to equations) +% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith +% equation contemporaneous restriction matrix where qi is the number of free parameters. +% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector +% of total original parameters and bi is a vector of free parameters. When no +% restrictions are imposed, we have Ui = I. There must be at least one free +% parameter left for the ith equation. +% nvar: number of endogeous variables +% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation +%---------------- +% b: sum(n0)-by-1 vector of A0 free parameters +% +% Tao Zha, February 2000. Revised, August 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +n0=n0(:); +n0cum = [0; cumsum(n0)]; +b=zeros(n0cum(end),1); +for kj = 1:nvar + b(n0cum(kj)+1:n0cum(kj+1))=Ui{kj}'*A0(:,kj); +end diff --git a/matlab/ms-sbvar/cstz/fn_tran_b2a.m b/matlab/ms-sbvar/cstz/fn_tran_b2a.m index bb06320388227951dc76ead549249679be3aae76..501b09326ad93c4d3702d106616ec790bb5b5fe9 100644 --- a/matlab/ms-sbvar/cstz/fn_tran_b2a.m +++ b/matlab/ms-sbvar/cstz/fn_tran_b2a.m @@ -1,41 +1,41 @@ -function A0 = fn_tran_b2a(b,Ui,nvar,n0) -% A0 = fn_tran_b2a(b,Ui,nvar,n0) -% Transform free parameters b's to A0. Note: columns correspond to equations -% -% b: sum(n0)-by-1 vector of A0 free parameters -% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith -% equation contemporaneous restriction matrix where qi is the number of free parameters. -% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector -% of total original parameters and bi is a vector of free parameters. When no -% restrictions are imposed, we have Ui = I. There must be at least one free -% parameter left for the ith equation. -% nvar: number of endogeous variables -% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation -%---------------- -% A0: nvar-by-nvar, contempareous matrix (columns correspond to equations) -% -% Tao Zha, February 2000. Revised, August 2000. - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -b=b(:); n0=n0(:); -A0 = zeros(nvar); -n0cum = [0; cumsum(n0)]; -for kj = 1:nvar - A0(:,kj) = Ui{kj}*b(n0cum(kj)+1:n0cum(kj+1)); -end +function A0 = fn_tran_b2a(b,Ui,nvar,n0) +% A0 = fn_tran_b2a(b,Ui,nvar,n0) +% Transform free parameters b's to A0. Note: columns correspond to equations +% +% b: sum(n0)-by-1 vector of A0 free parameters +% Ui: nvar-by-1 cell. In each cell, nvar-by-qi orthonormal basis for the null of the ith +% equation contemporaneous restriction matrix where qi is the number of free parameters. +% With this transformation, we have ai = Ui*bi or Ui'*ai = bi where ai is a vector +% of total original parameters and bi is a vector of free parameters. When no +% restrictions are imposed, we have Ui = I. There must be at least one free +% parameter left for the ith equation. +% nvar: number of endogeous variables +% n0: nvar-by-1, ith element represents the number of free A0 parameters in ith equation +%---------------- +% A0: nvar-by-nvar, contempareous matrix (columns correspond to equations) +% +% Tao Zha, February 2000. Revised, August 2000. + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +b=b(:); n0=n0(:); +A0 = zeros(nvar); +n0cum = [0; cumsum(n0)]; +for kj = 1:nvar + A0(:,kj) = Ui{kj}*b(n0cum(kj)+1:n0cum(kj+1)); +end diff --git a/matlab/ms-sbvar/cstz/fn_varoots.m b/matlab/ms-sbvar/cstz/fn_varoots.m index 0b7ae35c601e625cc0d619f58f81aa324c056545..9020b216844e133d2a38c50ad94f75cb57592763 100644 --- a/matlab/ms-sbvar/cstz/fn_varoots.m +++ b/matlab/ms-sbvar/cstz/fn_varoots.m @@ -1,46 +1,46 @@ -function rootsinv = fn_varoots(Bhat,nvar,lags) -% -% Using eigenvalues to find the inverse of all roots associated with the VAR proceess: -% y_t' = C + y_{t-1}'*B_1 + ... + Y_{t-p}'*B_p + u_t'. -% where columns correspond to equations. See also Judge (1), pp.753-755 where rows correspond to equations. -% Bhat: ncoef-by-nvar where ncoef=nvar*lags+nexo and nvar is the number of endogenous variables. -% Columns corresponds to equations with -% ncoef=[nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term] -% ..., nvar coef in the last lag, and nexo coefficients. -% Note that entries in the rows of Bhat that > nvar*lags are irrelevant. -% nvar: number of endogenous variables. -% lags: number of lags. -%------- -% rootsinv: a vector of nvar*lags inverse roots. When > 1, explosive. When all < 1, stationary. -% -% Tao Zha, September 2000 - -% Copyright (C) 2000-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -if size(Bhat,1)<nvar*lags - disp(' ') - warning('Make sure that Bhat has at least nvar*lags rows') - return -end - -%--------- Strack the VAR(p) to the VAR(1) with z_t = Az_{t-1}. -% -A1 = diag(ones(nvar*(lags-1),1)); -A2 = [A1 zeros(nvar*(lags-1),nvar)]; -A = [Bhat(1:nvar*lags,:)'; A2]; -rootsinv=eig(A); +function rootsinv = fn_varoots(Bhat,nvar,lags) +% +% Using eigenvalues to find the inverse of all roots associated with the VAR proceess: +% y_t' = C + y_{t-1}'*B_1 + ... + Y_{t-p}'*B_p + u_t'. +% where columns correspond to equations. See also Judge (1), pp.753-755 where rows correspond to equations. +% Bhat: ncoef-by-nvar where ncoef=nvar*lags+nexo and nvar is the number of endogenous variables. +% Columns corresponds to equations with +% ncoef=[nvar for 1st lag, ..., nvar for last lag, other exogenous terms, const term] +% ..., nvar coef in the last lag, and nexo coefficients. +% Note that entries in the rows of Bhat that > nvar*lags are irrelevant. +% nvar: number of endogenous variables. +% lags: number of lags. +%------- +% rootsinv: a vector of nvar*lags inverse roots. When > 1, explosive. When all < 1, stationary. +% +% Tao Zha, September 2000 + +% Copyright (C) 2000-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +if size(Bhat,1)<nvar*lags + disp(' ') + warning('Make sure that Bhat has at least nvar*lags rows') + return +end + +%--------- Strack the VAR(p) to the VAR(1) with z_t = Az_{t-1}. +% +A1 = diag(ones(nvar*(lags-1),1)); +A2 = [A1 zeros(nvar*(lags-1),nvar)]; +A = [Bhat(1:nvar*lags,:)'; A2]; +rootsinv=eig(A); diff --git a/matlab/ms-sbvar/cstz/sye.m b/matlab/ms-sbvar/cstz/sye.m index 60ddd60962381bec4f323b428070098eabd4e036..b4a24326e4a6808800aeb8c6d34345ae9dfee180 100644 --- a/matlab/ms-sbvar/cstz/sye.m +++ b/matlab/ms-sbvar/cstz/sye.m @@ -1,100 +1,100 @@ -function [Bh,e,xtx,xty,phi,y,ncoe,xr] = sye(z,lags) -% Now [Bh,e,xtx,xty,phi,y,ncoe,xr] = sye(z,lags) -% Old: [Bh,e,xtx,xty,phi,y,ncoe,Sigu,xtxinv] = sye(z,lags) -% -% Estimate a system of equations in the form of y(T*nvar) = XB + u, -% X--phi: T*k, B: k*nvar; where T=sp-lags, k=ncoe, -% -% z: (T+lags)-by-(nvar+1) raw data matrix (nvar of variables + constant). -% lags: number of lags -%-------------------- -% Bh: k-by-nvar estimated reduced-form parameter; column: nvar; -% row: k=ncoe=[nvar for 1st lag, ..., nvar for last lag, const] -% e: estimated residual e = y -xBh, T-by-nvar -% xtx: X'X: k-by-k -% xty: X'Y: k-by-nvar -% phi: X; T-by-k; column: [nvar for 1st lag, ..., nvar for last lag, const] -% y: Y: T-by-nvar -% ncoe: number of coeffcients per equation: nvar*lags + 1 -% xr: the economy size (k-by-k) in qr(phi) so that xr=chol(X'*X) -% Sigu: e'*e: nvar-by-nvar. Note, not divided (undivided) by "nobs" -% xtxinv: inv(X'X): k-by-k -% -% See also syed.m (allowing for more predetermined terms) which has not -% been yet updated as "sye.m". -% -% Note, "lags" is something I changed recently, so it may not be compatible -% with old programs, 10/15/98 by TAZ. -% -% Revised, 5/2/99. Replaced outputs Sigu and xtxinv with xr so that previous -% programs may be incompatible. - -% Copyright (C) 1998-2011 Tao Zha -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see <http://www.gnu.org/licenses/>. - -% ** setup of orders and lengths ** -[sp,nvar] = size(z); % sp: sample period T include lags -nvar = nvar-1; % -1: takes out the counting of constant - -ess = sp-lags; % effective sample size -sb = lags+1; % sample beginning -sl = sp; % sample last period -ncoe = nvar*lags + 1; % with constant - -% ** construct X for Y = X*B + U where phi = X ** -x = z(:,1:nvar); -C = z(:,nvar+1); -phi = zeros(ess,ncoe); -phi(:,ncoe) = C(1:ess); -for k=1:lags, phi(:,nvar*(k-1)+1:nvar*k) = x(sb-k:sl-k,:); end -% row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, const] -% Thus, # of columns is nvar*lags+1 = ncoef. -% ** estimate: B, XTX, residuals ** -y = x(sb:sl,:); -% -%**** The following, though stable, is too slow ***** -% [u d v]=svd(phi,0); %trial -% %xtx = phi'*phi; % X'X, k*k (ncoe*ncoe) -% vd=v.*(ones(size(v,2),1)*diag(d)'); %trial -% dinv = 1./diag(d); % inv(diag(d)) -% vdinv=v.*(ones(size(v,2),1)*dinv'); %trial -% xtx=vd*vd'; -% xtxinv = vdinv*vdinv'; -% %xty = phi'*y; % X'Y -% uy = u'*y; %trial -% xty = vd*uy; %trial -% %Bh = xtx\xty; %inv(X'X)*(X'Y), k*m (ncoe*nvar). -% Bh = xtxinv*xty; -% %e = y - phi*Bh; % from Y = XB + U, e: (T-lags)*nvar -% e = y - u*uy; -%**** The following, though stable, is too slow ***** - -%===== (Fast but perhaps less accurate) alternative to the above ========= -[xq,xr]=qr(phi,0); -xtx=xr'*xr; -xty=phi'*y; -Bh = xr\(xr'\xty); -e=y-phi*Bh; -%===== (Fast but perhaps less accurate) alternative to the above ========= - - -%* not numerically stable way of computing e'*e -%Sigu = y'*y-xty'*Bh; -%Sigu = y'*(eye(ess)-phi*xtxinv*phi')*y; % probablly better way, commented out - % by TZ, 2/28/98. See following -%Sigu = y'*(eye(ess)-u*u')*y; % I think this is the best, TZ, 2/28/98 - % Note, u*u'=x*inv(x'x)*x. +function [Bh,e,xtx,xty,phi,y,ncoe,xr] = sye(z,lags) +% Now [Bh,e,xtx,xty,phi,y,ncoe,xr] = sye(z,lags) +% Old: [Bh,e,xtx,xty,phi,y,ncoe,Sigu,xtxinv] = sye(z,lags) +% +% Estimate a system of equations in the form of y(T*nvar) = XB + u, +% X--phi: T*k, B: k*nvar; where T=sp-lags, k=ncoe, +% +% z: (T+lags)-by-(nvar+1) raw data matrix (nvar of variables + constant). +% lags: number of lags +%-------------------- +% Bh: k-by-nvar estimated reduced-form parameter; column: nvar; +% row: k=ncoe=[nvar for 1st lag, ..., nvar for last lag, const] +% e: estimated residual e = y -xBh, T-by-nvar +% xtx: X'X: k-by-k +% xty: X'Y: k-by-nvar +% phi: X; T-by-k; column: [nvar for 1st lag, ..., nvar for last lag, const] +% y: Y: T-by-nvar +% ncoe: number of coeffcients per equation: nvar*lags + 1 +% xr: the economy size (k-by-k) in qr(phi) so that xr=chol(X'*X) +% Sigu: e'*e: nvar-by-nvar. Note, not divided (undivided) by "nobs" +% xtxinv: inv(X'X): k-by-k +% +% See also syed.m (allowing for more predetermined terms) which has not +% been yet updated as "sye.m". +% +% Note, "lags" is something I changed recently, so it may not be compatible +% with old programs, 10/15/98 by TAZ. +% +% Revised, 5/2/99. Replaced outputs Sigu and xtxinv with xr so that previous +% programs may be incompatible. + +% Copyright (C) 1998-2011 Tao Zha +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +% ** setup of orders and lengths ** +[sp,nvar] = size(z); % sp: sample period T include lags +nvar = nvar-1; % -1: takes out the counting of constant + +ess = sp-lags; % effective sample size +sb = lags+1; % sample beginning +sl = sp; % sample last period +ncoe = nvar*lags + 1; % with constant + +% ** construct X for Y = X*B + U where phi = X ** +x = z(:,1:nvar); +C = z(:,nvar+1); +phi = zeros(ess,ncoe); +phi(:,ncoe) = C(1:ess); +for k=1:lags, phi(:,nvar*(k-1)+1:nvar*k) = x(sb-k:sl-k,:); end +% row: T-lags; column: [nvar for 1st lag, ..., nvar for last lag, const] +% Thus, # of columns is nvar*lags+1 = ncoef. +% ** estimate: B, XTX, residuals ** +y = x(sb:sl,:); +% +%**** The following, though stable, is too slow ***** +% [u d v]=svd(phi,0); %trial +% %xtx = phi'*phi; % X'X, k*k (ncoe*ncoe) +% vd=v.*(ones(size(v,2),1)*diag(d)'); %trial +% dinv = 1./diag(d); % inv(diag(d)) +% vdinv=v.*(ones(size(v,2),1)*dinv'); %trial +% xtx=vd*vd'; +% xtxinv = vdinv*vdinv'; +% %xty = phi'*y; % X'Y +% uy = u'*y; %trial +% xty = vd*uy; %trial +% %Bh = xtx\xty; %inv(X'X)*(X'Y), k*m (ncoe*nvar). +% Bh = xtxinv*xty; +% %e = y - phi*Bh; % from Y = XB + U, e: (T-lags)*nvar +% e = y - u*uy; +%**** The following, though stable, is too slow ***** + +%===== (Fast but perhaps less accurate) alternative to the above ========= +[xq,xr]=qr(phi,0); +xtx=xr'*xr; +xty=phi'*y; +Bh = xr\(xr'\xty); +e=y-phi*Bh; +%===== (Fast but perhaps less accurate) alternative to the above ========= + + +%* not numerically stable way of computing e'*e +%Sigu = y'*y-xty'*Bh; +%Sigu = y'*(eye(ess)-phi*xtxinv*phi')*y; % probablly better way, commented out + % by TZ, 2/28/98. See following +%Sigu = y'*(eye(ess)-u*u')*y; % I think this is the best, TZ, 2/28/98 + % Note, u*u'=x*inv(x'x)*x.