Description of bvar

0001 function [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags)
0002 %function [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags)
0003 % bvar_toolbox  Routines shared between BVAR methods
0004 % Computes several things for the estimations of a BVAR(nlags)
0005 %
0006 % INPUTS:
0007 %     nlags: number of lags
0008 %
0009 % OUTPUTS:
0010 %    ny:            number of endogenous variables
0011 %    nx:            number of exogenous variables (equal to zero, or one if a
0012 %                   constant term is included)
0013 %    posterior:     a structure describing the posterior distribution (which is
0014 %                   normal-Inverse-Wishart)
0015 %                   Its fields are:
0016 %                   - df: degrees of freedom of the inverse-Wishart distribution
0017 %                   - S: matrix parameter for the inverse-Wishart distribution
0018 %                   - XXi: first component of the VCV of the matrix-normal
0019 %                     distribution (the other one being drawn from the
0020 %                     inverse-Wishart)
0021 %                   - PhiHat: mean of the matrix-normal distribution
0022 %    prior:         a structure describing the prior distribution
0023 %                   Its fields are the same than for the posterior
0024 %    forecast_data: a structure containing data useful for forecasting
0025 %                   Its fields are:
0026 %                   - initval: a nlags*ny matrix containing the "nlags" last
0027 %                     observations of the sample (i.e. before options_.nobs)
0028 %                   - xdata: a matrix containing the future exogenous for
0029 %                     forecasting, of size options_.forecast*nx (actually only
0030 %                     contains "1" values for the constant term if nx ~= 0)
0031 %                   - realized_val: only non-empty if options_.nobs doesn't point
0032 %                     to the end of sample
0033 %                     In that case, contains values of endogenous variables after
0034 %                     options_.nobs and up to the end of the sample
0035 %                   - realized_xdata: contains values of exogenous variables after
0036 %                     options_.nobs and up to the end of the sample (actually only
0037 %                     contains "1" values for the constant term if nx ~= 0)
0038 %
0039 % SPECIAL REQUIREMENTS:
0040 %    This function uses the following Dynare options:
0041 %    - datafile, first_obs, varobs, xls_sheet, xls_range, nobs, presample
0042 %    - bvar_prior_{tau,decay,lambda,mu,omega,flat,train}
0043 
0044 % Copyright (C) 2003-2007 Christopher Sims
0045 % Copyright (C) 2007-2012 Dynare Team
0046 %
0047 % This file is part of Dynare.
0048 %
0049 % Dynare is free software: you can redistribute it and/or modify
0050 % it under the terms of the GNU General Public License as published by
0051 % the Free Software Foundation, either version 3 of the License, or
0052 % (at your option) any later version.
0053 %
0054 % Dynare is distributed in the hope that it will be useful,
0055 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0056 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0057 % GNU General Public License for more details.
0058 %
0059 % You should have received a copy of the GNU General Public License
0060 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0061 
0062 global options_
0063 
0064 % Load dataset
0065 dataset = read_variables(options_.datafile, options_.varobs, [], options_.xls_sheet, options_.xls_range);
0066 options_ = set_default_option(options_, 'nobs', size(dataset,1)-options_.first_obs+1);
0067 
0068 % Parameters for prior
0069 if options_.first_obs + options_.presample <= nlags
0070     error('first_obs+presample should be > nlags (for initializing the VAR)')
0071 end
0072 
0073 train = options_.bvar_prior_train;
0074 
0075 if options_.first_obs + options_.presample - train <= nlags
0076     error('first_obs+presample-train should be > nlags (for initializating the VAR)')
0077 end
0078 
0079 idx = options_.first_obs+options_.presample-train-nlags:options_.first_obs+options_.nobs-1;
0080 
0081 % Prepare dataset
0082 if options_.loglinear && ~options_.logdata
0083     dataset = log(dataset);
0084 end
0085 if options_.prefilter
0086     dataset(idx,:) = dataset(idx,:) - ones(length(idx),1)*mean(dataset(idx,:));
0087 end
0088 
0089 mnprior.tight = options_.bvar_prior_tau;
0090 mnprior.decay = options_.bvar_prior_decay;
0091 
0092 % Use only initializations lags for the variance prior
0093 vprior.sig = std(dataset(options_.first_obs+options_.presample-nlags:options_.first_obs+options_.presample,:))';
0094 vprior.w = options_.bvar_prior_omega;
0095 
0096 lambda = options_.bvar_prior_lambda;
0097 mu = options_.bvar_prior_mu;
0098 flat = options_.bvar_prior_flat;
0099 
0100 ny = size(dataset, 2);
0101 if options_.prefilter || options_.noconstant
0102     nx = 0;
0103 else
0104     nx = 1;
0105 end
0106 
0107 [ydum, xdum, pbreaks] = varprior(ny, nx, nlags, mnprior, vprior);
0108 
0109 ydata = dataset(idx, :);
0110 T = size(ydata, 1);
0111 xdata = ones(T,nx);
0112 
0113 % Posterior density
0114 var = rfvar3([ydata; ydum], nlags, [xdata; xdum], [T; T+pbreaks], lambda, mu);
0115 Tu = size(var.u, 1);
0116 
0117 posterior.df = Tu - ny*nlags - nx - flat*(ny+1);
0118 posterior.S = var.u' * var.u;
0119 posterior.XXi = var.xxi;
0120 posterior.PhiHat = var.B;
0121 
0122 % Prior density
0123 Tp = train + nlags;
0124 if nx
0125     xdata = xdata(1:Tp, :);
0126 else
0127     xdata = [];
0128 end
0129 varp = rfvar3([ydata(1:Tp, :); ydum], nlags, [xdata; xdum], [Tp; Tp + pbreaks], lambda, mu);
0130 Tup = size(varp.u, 1);
0131 
0132 prior.df = Tup - ny*nlags - nx - flat*(ny+1);
0133 prior.S = varp.u' * varp.u;
0134 prior.XXi = varp.xxi;
0135 prior.PhiHat = varp.B;
0136 
0137 if prior.df < ny
0138     error('Too few degrees of freedom in the inverse-Wishart part of prior distribution. You should increase training sample size.')
0139 end
0140 
0141 % Add forecast informations
0142 if nargout >= 5
0143     forecast_data.xdata = ones(options_.forecast, nx);
0144     forecast_data.initval = ydata(end-nlags+1:end, :);
0145     if options_.first_obs + options_.nobs <= size(dataset, 1)
0146         forecast_data.realized_val = dataset(options_.first_obs+options_.nobs:end, :);
0147         forecast_data.realized_xdata = ones(size(forecast_data.realized_val, 1), nx);
0148     else
0149         forecast_data.realized_val = [];
0150     end
0151 end
0152 
0153 
0154 function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
0155 %function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
0156 % ydum, xdum:   dummy observation data that implement the prior
0157 % breaks:       vector of points in the dummy data after which new dummy obs's start
0158 %                   Set breaks=T+[0;breaks], ydata=[ydata;ydum], xdum=[xdata;xdum], where
0159 %                   actual data matrix has T rows, in preparing input for rfvar3
0160 % nv,nx,lags: VAR dimensions
0161 % mnprior.tight:Overall tightness of Minnesota prior
0162 % mnprior.decay:Standard deviations of lags shrink as lag^(-decay)
0163 % vprior.sig:   Vector of prior modes for diagonal elements of r.f. covariance matrix
0164 % vprior.w:     Weight on prior on vcv.  1 corresponds to "one dummy observation" weight
0165 %                   Should be an integer, and will be rounded if not.  vprior.sig is needed
0166 %                   to scale the Minnesota prior, even if the prior on sigma is not used itself.
0167 %                   Set vprior.w=0 to achieve this.
0168 % Note:         The original Minnesota prior treats own lags asymmetrically, and therefore
0169 %                   cannot be implemented entirely with dummy observations.  It is also usually
0170 %                   taken to include the sum-of-coefficients and co-persistence components
0171 %                   that are implemented directly in rfvar3.m.  The diagonal prior on v, combined
0172 %                   with sum-of-coefficients and co-persistence components and with the unit own-first-lag
0173 %                   prior mean generates larger prior variances for own than for cross-effects even in
0174 %                   this formulation, but here there is no way to shrink toward a set of unconstrained
0175 %                   univariate AR's.
0176 
0177 % Original file downloaded from:
0178 % http://sims.princeton.edu/yftp/VARtools/matlab/varprior.m
0179 
0180 if ~isempty(mnprior)
0181     xdum = zeros(lags+1,nx,lags,nv);
0182     ydum = zeros(lags+1,nv,lags,nv);
0183     for il = 1:lags
0184         ydum(il+1,:,il,:) = il^mnprior.decay*diag(vprior.sig);
0185     end
0186     ydum(1,:,1,:) = diag(vprior.sig);
0187     ydum = mnprior.tight*reshape(ydum,[lags+1,nv,lags*nv]);
0188     ydum = flipdim(ydum,1);
0189     xdum = mnprior.tight*reshape(xdum,[lags+1,nx,lags*nv]);
0190     xdum = flipdim(xdum,1);
0191     breaks = (lags+1)*[1:(nv*lags)]';
0192     lbreak = breaks(end);
0193 else
0194     ydum = [];
0195     xdum = [];
0196     breaks = [];
0197     lbreak = 0;
0198 end
0199 if ~isempty(vprior) && vprior.w>0
0200     ydum2 = zeros(lags+1,nv,nv);
0201     xdum2 = zeros(lags+1,nx,nv);
0202     ydum2(end,:,:) = diag(vprior.sig);
0203     for i = 1:vprior.w
0204         ydum = cat(3,ydum,ydum2);
0205         xdum = cat(3,xdum,xdum2);
0206         breaks = [breaks;(lags+1)*[1:nv]'+lbreak];
0207         lbreak = breaks(end);
0208     end
0209 end
0210 dimy = size(ydum);
0211 ydum = reshape(permute(ydum,[1 3 2]),dimy(1)*dimy(3),nv);
0212 xdum = reshape(permute(xdum,[1 3 2]),dimy(1)*dimy(3),nx);
0213 breaks = breaks(1:(end-1));
0214 
0215 
0216 function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
0217 %function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
0218 % This algorithm goes for accuracy without worrying about memory requirements.
0219 % ydata:   dependent variable data matrix
0220 % xdata:   exogenous variable data matrix
0221 % lags:    number of lags
0222 % breaks:  rows in ydata and xdata after which there is a break.  This allows for
0223 %          discontinuities in the data (e.g. war years) and for the possibility of
0224 %          adding dummy observations to implement a prior.  This must be a column vector.
0225 %          Note that a single dummy observation becomes lags+1 rows of the data matrix,
0226 %          with a break separating it from the rest of the data.  The function treats the
0227 %          first lags observations at the top and after each "break" in ydata and xdata as
0228 %          initial conditions.
0229 % lambda:  weight on "co-persistence" prior dummy observations.  This expresses
0230 %          belief that when data on *all* y's are stable at their initial levels, they will
0231 %          tend to persist at that level.  lambda=5 is a reasonable first try.  With lambda<0,
0232 %          constant term is not included in the dummy observation, so that stationary models
0233 %          with means equal to initial ybar do not fit the prior mean.  With lambda>0, the prior
0234 %          implies that large constants are unlikely if unit roots are present.
0235 % mu:      weight on "own persistence" prior dummy observation.  Expresses belief
0236 %          that when y_i has been stable at its initial level, it will tend to persist
0237 %          at that level, regardless of the values of other variables.  There is
0238 %          one of these for each variable.  A reasonable first guess is mu=2.
0239 %      The program assumes that the first lags rows of ydata and xdata are real data, not dummies.
0240 %      Dummy observations should go at the end, if any.  If pre-sample x's are not available,
0241 %      repeating the initial xdata(lags+1,:) row or copying xdata(lags+1:2*lags,:) into
0242 %      xdata(1:lags,:) are reasonable subsititutes.  These values are used in forming the
0243 %      persistence priors.
0244 
0245 % Original file downloaded from:
0246 % http://sims.princeton.edu/yftp/VARtools/matlab/rfvar3.m
0247 
0248 [T,nvar] = size(ydata);
0249 nox = isempty(xdata);
0250 if ~nox
0251     [T2,nx] = size(xdata);
0252 else
0253     T2 = T;
0254     nx = 0;
0255     xdata = zeros(T2,0);
0256 end
0257 % note that x must be same length as y, even though first part of x will not be used.
0258 % This is so that the lags parameter can be changed without reshaping the xdata matrix.
0259 if T2 ~= T, error('Mismatch of x and y data lengths'),end
0260 if nargin < 4
0261     nbreaks = 0;
0262     breaks = [];
0263 else
0264     nbreaks = length(breaks);
0265 end
0266 breaks = [0;breaks;T];
0267 smpl = [];
0268 for nb = 1:nbreaks+1
0269     smpl = [smpl;[breaks(nb)+lags+1:breaks(nb+1)]'];
0270 end
0271 Tsmpl = size(smpl,1);
0272 X = zeros(Tsmpl,nvar,lags);
0273 for is = 1:length(smpl)
0274     X(is,:,:) = ydata(smpl(is)-(1:lags),:)';
0275 end
0276 X = [X(:,:) xdata(smpl,:)];
0277 y = ydata(smpl,:);
0278 % Everything now set up with input data for y=Xb+e
0279 
0280 % Add persistence dummies
0281 if lambda ~= 0 || mu > 0
0282     ybar = mean(ydata(1:lags,:),1);
0283     if ~nox
0284         xbar = mean(xdata(1:lags,:),1);
0285     else
0286         xbar = [];
0287     end
0288     if lambda ~= 0
0289         if lambda>0
0290             xdum = lambda*[repmat(ybar,1,lags) xbar];
0291         else
0292             lambda = -lambda;
0293             xdum = lambda*[repmat(ybar,1,lags) zeros(size(xbar))];
0294         end
0295         ydum = zeros(1,nvar);
0296         ydum(1,:) = lambda*ybar;
0297         y = [y;ydum];
0298         X = [X;xdum];
0299     end
0300     if mu>0
0301         xdum = [repmat(diag(ybar),1,lags) zeros(nvar,nx)]*mu;
0302         ydum = mu*diag(ybar);
0303         X = [X;xdum];
0304         y = [y;ydum];
0305     end
0306 end
0307 
0308 % Compute OLS regression and residuals
0309 [vl,d,vr] = svd(X,0);
0310 di = 1./diag(d);
0311 B = (vr.*repmat(di',nvar*lags+nx,1))*vl'*y;
0312 u = y-X*B;
0313 xxi = vr.*repmat(di',nvar*lags+nx,1);
0314 xxi = xxi*xxi';
0315 
0316 var.B = B;
0317 var.u = u;
0318 var.xxi = xxi;
bvar_toolbox

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
