v4: added BVAR code for computing marginal density and forecasting

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1331 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
sebastien 2007-06-25 22:42:30 +00:00
parent ec4d29935d
commit 86360d9eb4
15 changed files with 2774 additions and 2245 deletions

62
matlab/bvar_density.m Normal file
View File

@ -0,0 +1,62 @@
function bvar_density(maxnlags)
for nlags = 1:maxnlags
[ny, nx, posterior, prior] = bvar_toolbox(nlags);
posterior_int = matrictint(posterior.S, posterior.df, posterior.XXi);
prior_int = matrictint(prior.S, prior.df, prior.XXi);
lik_nobs = posterior.df - prior.df;
log_dnsty = posterior_int - prior_int - 0.5*ny*lik_nobs*log(2*pi);
disp(' ')
fprintf('The marginal log density of the BVAR(%g) model is equal to %10.4f\n', ...
nlags, log_dnsty);
disp(' ')
end
function w = matrictint(S, df, XXi)
% Computes the integral of the kernel of the PDF of a
% normal-inverse-Wishart distribution.
%
% S: parameter of inverse-Wishart distribution
% df: number of degrees of freedom of inverse-Wishart distribution
% XXi: first component of VCV matrix of matrix-normal distribution
%
% Computes the integral over (Phi, Sigma) of:
%
% det(Sigma)^(-k/2)*exp(-0.5*Tr((Phi-PhiHat)'*(XXi)^(-1)*(Phi-PhiHat)*Sigma^(-1)))*
% det(Sigma)^((df+ny+1)/2)*exp(-0.5*Tr(Sigma^(-1)*S))
%
% (where k is the dimension of XXi and ny is the dimension of S and
% Sigma)
k=size(XXi,1);
ny=size(S,1);
[cx,p]=chol(XXi);
[cs,q]=chol(S);
if any(diag(cx)<100*eps)
error('singular XXi')
end
if any(diag(cs<100*eps))
error('singular S')
end
% Matrix-normal component
w1 = 0.5*k*ny*log(2*pi)+ny*sum(log(diag(cx)));
% Inverse-Wishart component
w2 = -df*sum(log(diag(cs))) + 0.5*df*ny*log(2) + ny*(ny-1)*0.25*log(pi) + ggammaln(ny, df);
w = w1 + w2;
function lgg = ggammaln(m, df)
if df <= (m-1)
error('too few df in ggammaln')
else
garg = 0.5*(df+(0:-1:1-m));
lgg = sum(gammaln(garg));
end

152
matlab/bvar_forecast.m Normal file
View File

@ -0,0 +1,152 @@
function bvar_forecast(nlags)
global options_
options_ = set_default_option(options_, 'bvar_replic', 2000);
if options_.forecast == 0
error('bvar_forecast: you must specify "forecast" option')
end
[ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags);
sims_no_shock = NaN(options_.forecast, ny, options_.bvar_replic);
sims_with_shocks = NaN(options_.forecast, ny, options_.bvar_replic);
S_inv_chol = chol(inv(posterior.S));
XXi_lower_chol = chol(posterior.XXi, 'lower');
k = ny*nlags+nx;
for d = 1:options_.bvar_replic
Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_chol);
Sigma_lower_chol = chol(Sigma, 'lower');
Phi = rand_matrix_normal(k, ny, posterior.PhiHat, XXi_lower_chol, Sigma_lower_chol);
% Without shocks
lags_data = forecast_data.initval;
for t = 1:options_.forecast
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
y = X * Phi;
lags_data = [ lags_data(2:end, :); y ];
sims_no_shock(t, :, d) = y;
end
% With shocks
lags_data = forecast_data.initval;
for t = 1:options_.forecast
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
shock = (Sigma_lower_chol * randn(ny, 1))';
y = X * Phi + shock;
lags_data = [ lags_data(2:end, :); y ];
sims_with_shocks(t, :, d) = y;
end
end
% Plot graphs
sims_no_shock_mean = mean(sims_no_shock, 3);
sort_idx = round((0.5 + [-options_.conf_sig, options_.conf_sig]/2) * options_.bvar_replic);
sims_no_shock_sort = sort(sims_no_shock, 3);
sims_no_shock_down_conf = sims_no_shock_sort(:, :, sort_idx(1));
sims_no_shock_up_conf = sims_no_shock_sort(:, :, sort_idx(2));
sims_with_shocks_sort = sort(sims_with_shocks, 3);
sims_with_shocks_down_conf = sims_with_shocks_sort(:, :, sort_idx(1));
sims_with_shocks_up_conf = sims_with_shocks_sort(:, :, sort_idx(2));
dynare_graph_init(sprintf('BVAR forecasts (nlags = %d)', nlags), ny, {'b-' 'g-' 'g-' 'r-' 'r-'});
for i = 1:ny
dynare_graph([ sims_no_shock_mean(:, i) ...
sims_no_shock_up_conf(:, i) sims_no_shock_down_conf(:, i) ...
sims_with_shocks_up_conf(:, i) sims_with_shocks_down_conf(:, i) ], ...
options_.varobs(i, :));
end
dynare_graph_close;
% Compute RMSE
if ~isempty(forecast_data.realized_val)
sq_err_cumul = zeros(1, ny);
lags_data = forecast_data.initval;
for t = 1:size(forecast_data.realized_val, 1)
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.realized_xdata(t, :) ];
y = X * Phi;
lags_data = [ lags_data(2:end, :); y ];
sq_err_cumul = sq_err_cumul + (y - forecast_data.realized_val(t, :)) .^ 2;
end
rmse = sqrt(sq_err_cumul / size(sq_err_cumul, 1));
fprintf('RMSE of BVAR(%d):\n', nlags);
for i = 1:size(options_.varobs, 1)
fprintf('%s: %10.4f\n', options_.varobs(i, :), rmse(i));
end
end
function G = rand_inverse_wishart(m, v, H_inv_chol)
% rand_inverse_wishart Pseudo random matrices drawn from an
% inverse Wishart distribution
%
% G = rand_inverse_wishart(m, v, H_inv_chol)
%
% Returns an m-by-m matrix drawn from an inverse-Wishart distribution.
%
% m: dimension of G and H_inv_chol.
% v: degrees of freedom, greater or equal than m.
% H_inv_chol: (upper) cholesky decomposition of the inverse of the
% matrix parameter.
% The (upper) cholesky of the inverse is requested here
% in order to avoid to recompute it at every random draw.
% H_inv_chol = chol(inv(H))
%
% In other words:
% G ~ IW(m, v, H) where H = inv(H_inv_chol'*H_inv_chol)
% or, equivalently, using the correspondence between Wishart and
% inverse-Wishart:
% inv(G) ~ W(m, v, S) where S = H_inv_chol'*H_inv_chol = inv(H)
X = NaN(v, m);
for i = 1:v
X(i, :) = randn(1, m) * H_inv_chol;
end
% At this point, X'*X is Wishart distributed
% G = inv(X'*X);
% Rather compute inv(X'*X) using the SVD
[U,S,V] = svd(X, 0);
SSi = 1 ./ (diag(S) .^ 2);
G = (V .* repmat(SSi', m, 1)) * V';
function B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol)
% rand_matrix_normal Pseudo random matrices drawn from a
% matrix-normal distribution
%
% B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol)
%
% Returns an n-by-p matrix drawn from a Matrix-normal distribution
% Same notations than: http://en.wikipedia.org/wiki/Matrix_normal_distribution
% M is the mean, n-by-p matrix
% Omega_lower_chol is n-by-n, lower Cholesky decomposition of Omega
% (Omega_lower_chol = chol(Omega, 'lower'))
% Sigma_lower_chol is p-by-p, lower Cholesky decomposition of Sigma
% (Sigma_lower_chol = chol(Sigma, 'lower'))
%
% Equivalent to vec(B) ~ N(vec(Mu), kron(Sigma, Omega))
%
B1 = randn(n * p, 1);
B2 = kron(Sigma_lower_chol, Omega_lower_chol) * B1;
B3 = reshape(B2, n, p);
B = B3 + M;

277
matlab/bvar_toolbox.m Normal file
View File

@ -0,0 +1,277 @@
function [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags)
% bvar_toolbox Routines shared between BVAR methods
%
% [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags)
%
% Computes several things for the estimations of a BVAR(nlags):
% ny: number of endogenous variables
% nx: number of exogenous variables (actually equal to one, that is
% the constant term)
% posterior: a structure describing the posterior distribution (which is
% normal-Inverse-Wishart)
% Its fields are:
% - df: degrees of freedom of the inverse-Wishart distribution
% - S: matrix parameter for the inverse-Wishart distribution
% - XXi: first component of the VCV of the matrix-normal
% distribution (the other one being drawn from the
% inverse-Wishart)
% - PhiHat: mean of the matrix-normal distribution
% prior: a structure describing the prior distribution
% Its fields are the same than for the posterior
% forecast: a structure containing data useful for forecasting
% Its fields are:
% - initval: a nlags*ny vector containing the "nlags" last
% observations of the sample (i.e. before options_.nobs)
% - xdata: a vector containing the future exogenous for
% forecasting, of size options_.forecast*nx (actually only
% contains "1" values for the constant term)
% - realized_val: only non-empty if options_.nobs doesn't point
% to the end of sample
% In that case, contains values of endogenous variables after
% options_.nobs and up to the end of the sample
% - realized_xdata: contains values of exogenous variables after
% options_.nobs and up to the end of the sample (actually only
% contains "1" values)
%
% This function uses the following Dynare options:
% - datafile, first_obs, varobs, xls_sheet, xls_range, nobs, presample
% - bvar_prior_{tau,decay,lambda,mu,omega,flat,train}
global options_
% Prepare dataset
dataset = read_variables(options_.datafile, options_.varobs, [], options_.xls_sheet, options_.xls_range);
options_ = set_default_option(options_, 'nobs', size(dataset,1)-options_.first_obs+1);
% Parameters for prior
options_ = set_default_option(options_, 'bvar_prior_tau', 3);
options_ = set_default_option(options_, 'bvar_prior_decay', 0.5);
options_ = set_default_option(options_, 'bvar_prior_lambda', 5);
options_ = set_default_option(options_, 'bvar_prior_mu', 2);
options_ = set_default_option(options_, 'bvar_prior_omega', 1);
options_ = set_default_option(options_, 'bvar_prior_flat', 0);
options_ = set_default_option(options_, 'bvar_prior_train', 0);
if options_.first_obs + options_.presample <= nlags
error('first_obs+presample should be > nlags (for initializing the VAR)')
end
train = options_.bvar_prior_train;
if options_.first_obs + options_.presample - train <= nlags
error('first_obs+presample-train should be > nlags (for initializating the VAR)')
end
mnprior.tight = options_.bvar_prior_tau;
mnprior.decay = options_.bvar_prior_decay;
% Use only initializations lags for the variance prior
vprior.sig = std(dataset(options_.first_obs+options_.presample-nlags:options_.first_obs+options_.presample-1,:))';
vprior.w = options_.bvar_prior_omega;
lambda = options_.bvar_prior_lambda;
mu = options_.bvar_prior_mu;
flat = options_.bvar_prior_flat;
ny = size(dataset, 2);
nx = 1;
[ydum, xdum, pbreaks] = varprior(ny, nx, nlags, mnprior, vprior);
ydata = dataset(options_.first_obs+options_.presample-train-nlags:options_.first_obs+options_.nobs-1, :);
T = size(ydata, 1);
xdata = ones(T, 1);
% Posterior density
var = rfvar3([ydata; ydum], nlags, [xdata; xdum], [T; T+pbreaks], lambda, mu);
Tu = size(var.u, 1);
posterior.df = Tu - ny*nlags - nx - flat*(ny+1);
posterior.S = var.u' * var.u;
posterior.XXi = var.xxi;
posterior.PhiHat = var.B;
% Prior density
Tp = train + nlags;
varp = rfvar3([ydata(1:Tp, :); ydum], nlags, [xdata(1:Tp, :); xdum], ...
[Tp; Tp + pbreaks], lambda, mu);
Tup = size(varp.u, 1);
prior.df = Tup - ny*nlags - nx - flat*(ny+1);
prior.S = varp.u' * varp.u;
prior.XXi = varp.xxi;
prior.PhiHat = varp.B;
if prior.df < ny
error('Too few degrees of freedom in the inverse-Wishart part of prior distribution. You should increase training sample size.')
end
% Add forecast informations
if nargout >= 5
forecast_data.xdata = ones(options_.forecast, nx);
forecast_data.initval = ydata(end-nlags+1:end, :);
if options_.first_obs + options_.nobs <= size(dataset, 1)
forecast_data.realized_val = dataset(options_.first_obs+options_.nobs:end, :);
forecast_data.realized_xdata = ones(size(forecast_data.realized_val, 1), 1);
else
forecast_data.realized_val = [];
end
end
function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
%function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
% ydum, xdum: dummy observation data that implement the prior
% breaks: vector of points in the dummy data after which new dummy obs's start
% Set breaks=T+[0;breaks], ydata=[ydata;ydum], xdum=[xdata;xdum], where
% actual data matrix has T rows, in preparing input for rfvar3
% nv,nx,lags: VAR dimensions
% mnprior.tight:Overall tightness of Minnesota prior
% mnprior.decay:Standard deviations of lags shrink as lag^(-decay)
% vprior.sig: Vector of prior modes for diagonal elements of r.f. covariance matrix
% vprior.w: Weight on prior on vcv. 1 corresponds to "one dummy observation" weight
% Should be an integer, and will be rounded if not. vprior.sig is needed
% to scale the Minnesota prior, even if the prior on sigma is not used itself.
% Set vprior.w=0 to achieve this.
% Note: The original Minnesota prior treats own lags asymmetrically, and therefore
% cannot be implemented entirely with dummy observations. It is also usually
% taken to include the sum-of-coefficients and co-persistence components
% that are implemented directly in rfvar3.m. The diagonal prior on v, combined
% with sum-of-coefficients and co-persistence components and with the unit own-first-lag
% prior mean generates larger prior variances for own than for cross-effects even in
% this formulation, but here there is no way to shrink toward a set of unconstrained
% univariate AR's.
% Author: C. Sims
if ~isempty(mnprior)
xdum = zeros(lags+1,nx,lags,nv);
ydum = zeros(lags+1,nv,lags,nv);
for il = 1:lags
ydum(il+1,:,il,:) = il^mnprior.decay*diag(vprior.sig);
end
ydum(1,:,1,:) = diag(vprior.sig);
ydum = mnprior.tight*reshape(ydum,[lags+1,nv,lags*nv]);
ydum = flipdim(ydum,1);
xdum = mnprior.tight*reshape(xdum,[lags+1,nx,lags*nv]);
xdum = flipdim(xdum,1);
breaks = (lags+1)*[1:(nv*lags)]';
lbreak = breaks(end);
else
ydum = [];
xdum = [];
breaks = [];
lbreak = 0;
end
if ~isempty(vprior) & vprior.w>0
ydum2 = zeros(lags+1,nv,nv);
xdum2 = zeros(lags+1,nx,nv);
ydum2(end,:,:) = diag(vprior.sig);
ydum = cat(3,ydum,ydum2);
xdum = cat(3,xdum,xdum2);
dimy = size(ydum);
ydum = reshape(permute(ydum,[1 3 2]),dimy(1)*dimy(3),nv);
xdum = reshape(permute(xdum,[1 3 2]),dimy(1)*dimy(3),nx);
breaks = [breaks;(lags+1)*[1:nv-1]'+lbreak];
end
function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
%function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
% This algorithm goes for accuracy without worrying about memory requirements.
% ydata: dependent variable data matrix
% xdata: exogenous variable data matrix
% lags: number of lags
% breaks: rows in ydata and xdata after which there is a break. This allows for
% discontinuities in the data (e.g. war years) and for the possibility of
% adding dummy observations to implement a prior. This must be a column vector.
% Note that a single dummy observation becomes lags+1 rows of the data matrix,
% with a break separating it from the rest of the data. The function treats the
% first lags observations at the top and after each "break" in ydata and xdata as
% initial conditions.
% lambda: weight on "co-persistence" prior dummy observations. This expresses
% belief that when data on *all* y's are stable at their initial levels, they will
% tend to persist at that level. lambda=5 is a reasonable first try. With lambda<0,
% constant term is not included in the dummy observation, so that stationary models
% with means equal to initial ybar do not fit the prior mean. With lambda>0, the prior
% implies that large constants are unlikely if unit roots are present.
% mu: weight on "own persistence" prior dummy observation. Expresses belief
% that when y_i has been stable at its initial level, it will tend to persist
% at that level, regardless of the values of other variables. There is
% one of these for each variable. A reasonable first guess is mu=2.
% The program assumes that the first lags rows of ydata and xdata are real data, not dummies.
% Dummy observations should go at the end, if any. If pre-sample x's are not available,
% repeating the initial xdata(lags+1,:) row or copying xdata(lags+1:2*lags,:) into
% xdata(1:lags,:) are reasonable subsititutes. These values are used in forming the
% persistence priors.
% Code written by Christopher Sims. This version 6/15/03.
[T,nvar] = size(ydata);
nox = isempty(xdata);
if ~nox
[T2,nx] = size(xdata);
else
T2 = T;
nx = 0;
xdata = zeros(T2,0);
end
% note that x must be same length as y, even though first part of x will not be used.
% This is so that the lags parameter can be changed without reshaping the xdata matrix.
if T2 ~= T, error('Mismatch of x and y data lengths'),end
if nargin < 4
nbreaks = 0;
breaks = [];
else
nbreaks = length(breaks);
end
breaks = [0;breaks;T];
smpl = [];
for nb = 1:nbreaks+1
smpl = [smpl;[breaks(nb)+lags+1:breaks(nb+1)]'];
end
Tsmpl = size(smpl,1);
X = zeros(Tsmpl,nvar,lags);
for is = 1:length(smpl)
X(is,:,:) = ydata(smpl(is)-(1:lags),:)';
end
X = [X(:,:) xdata(smpl,:)];
y = ydata(smpl,:);
% Everything now set up with input data for y=Xb+e
% Add persistence dummies
if lambda ~= 0 | mu > 0
ybar = mean(ydata(1:lags,:),1);
if ~nox
xbar = mean(xdata(1:lags,:),1);
else
xbar = [];
end
if lambda ~= 0
if lambda>0
xdum = lambda*[repmat(ybar,1,lags) xbar];
else
lambda = -lambda;
xdum = lambda*[repmat(ybar,1,lags) zeros(size(xbar))];
end
ydum = zeros(1,nvar);
ydum(1,:) = lambda*ybar;
y = [y;ydum];
X = [X;xdum];
end
if mu>0
xdum = [repmat(diag(ybar),1,lags) zeros(nvar,nx)]*mu;
ydum = mu*diag(ybar);
X = [X;xdum];
y = [y;ydum];
end
end
% Compute OLS regression and residuals
[vl,d,vr] = svd(X,0);
di = 1./diag(d);
B = (vr.*repmat(di',nvar*lags+nx,1))*vl'*y;
u = y-X*B;
xxi = vr.*repmat(di',nvar*lags+nx,1);
xxi = xxi*xxi';
var.B = B;
var.u = u;
var.xxi = xxi;

View File

@ -1,35 +0,0 @@
function w=matrictint(S,XXi,T)
%function w=matrictint(S,XXi,T)
% S: usually sample cross product matrix of LS residuals
% XXi: inv(X'X) matrix for rhs variables
% T: number of observations
% w: log of integrated posterior for SUR or RF VAR with det(Sigma)^(-(m+1)/2) Jeffreys-like prior
% To get the log of the integral of the likelihood for a VAR with T observations,
% k rhs variables in each equation, and m equations, set T=T-m-1 and subtract .5*m*(m+1)*log(2*pi).
% We are integrating the exponential of -.5*T*m*log(2*pi)-.5*(T+m+1)*log(det(Sigma))-.5*trace(Sigma\S(beta)).
k=size(XXi,1);
m=size(S,1);
[cx,p]=chol(XXi);
[cs,q]=chol(S);
%cx=cschol(XXi);
%cs=cschol(S);
if any(diag(cx)<100*eps)
error('singular XXi')
end
if any(diag(cs<100*eps))
error('singular S')
end
w=(-T+k+(m-1)/2)*m*.5*log(pi)-(T-k)*sum(log(diag(cs)))+m*sum(log(diag(cx)))+ggammaln(m,(T-k)/2);
function lgg=ggammaln(m,ndf)
%function gg=ggamma(m,ndf)
% From 8.2.22 on p.427 of Box and Tiao, this is the log of generalized
% gamma divided by gamma(.5)^(.5*m*(m-1))
if ndf<=(m-1)/2
error('too few df in ggammaln')
else
%lgg=.5*m*(m-1)*gammaln(.5); % normalizing factor not used in Wishart integral
garg=ndf+.5*(0:-1:1-m);
lgg=sum(gammaln(garg));
end

View File

@ -1,50 +0,0 @@
function w=mgnldnsty(ydata,lags,xdata,breaks,lambda,mu,mnprior,vprior,train,flat)
%function w=mgnldnsty(ydata,lags,xdata,breaks,lambda,mu,mnprior,vprior,train,flat)
% ydata: endogenous variable data matrix, including initial condition dates.
% xdata: exogenous variable data matrix, including initial condition dates.
% breaks: breaks in the data. The first lags data points after a break are used
% as new initial conditions, not data points for the fit.
% lambda: weight on the co-persistence prior dummy observation. (5 is reasonable)
% lambda>0 => x variables included; lambda<0 => x variables excluded;
% mu: weight on variable-by-variable sum of coeffs dummy obs. (1 is reasonable)
% mnprior.tight:weight on the Minnesota prior dummies. Prior std dev on first lag is
% 1/mnprior.tight
% mnprior.decay:prior std dev on lag j is 1/j^decay
% vprior.sig: vector of nv prior std dev's of equation shocks
% vprior.w: weight on vcv dummies. (1 is reasonable; higher values tighten up.)
% train: If present and non-zero, this is the point in the sample at which the
% "training sample" ends. Prior x likelihood to this point is weighted to
% integrate to 1, and therefore is treated as if it were itself the prior.
% To do a pure training sample prior, set lambda=mu=0, mnprior=vprior=[],
% train>lags.
%
%flat: Even with lambda=mu=0, vprior=mnprior=[], det(Sigma)^(-(nv+1)/2) is used
% as a "prior", unless flat=1. flat, if present, must be 1 or 0.
% flat=1 is likely not to work unless train is reasonably large.
if nargin<10,flat=0;end
[T,nv]=size(ydata);
[Tx,nx]=size(xdata);
if Tx ~= T, error('ydata and xdata length mismatch'),end
[ydum,xdum,pbreaks]=varprior(nv,nx,lags,mnprior,vprior);
var=rfvar3([ydata;ydum],lags,[xdata;xdum],[breaks;T;T+pbreaks],lambda,mu);
Tu=size(var.u,1);
w=matrictint(var.u'*var.u,var.xxi,Tu-flat*(nv+1))-flat*.5*nv*(nv+1)*log(2*pi);
if nargin>8
if ~isempty(train) & train>0
if train <= lags
error('end of training sample <= # of lags')
end
Tp=train;
tbreaks=breaks(find(breaks<train));
else
Tp=lags;
tbreaks=[];
end
else
Tp=lags;
tbreaks=[];
end
varp=rfvar3([ydata(1:Tp,:);ydum],lags,[xdata(1:Tp);xdum],[tbreaks;Tp;Tp+pbreaks],lambda,mu);
Tup=size(varp.u,1);
wp=matrictint(varp.u'*varp.u,varp.xxi,Tup-flat*(nv+1)/2)-flat*.5*nv*(nv+1)*log(2*pi);
w=w-wp;

View File

@ -1,106 +0,0 @@
function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
%function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
% This algorithm goes for accuracy without worrying about memory requirements.
% ydata: dependent variable data matrix
% xdata: exogenous variable data matrix
% lags: number of lags
% breaks: rows in ydata and xdata after which there is a break. This allows for
% discontinuities in the data (e.g. war years) and for the possibility of
% adding dummy observations to implement a prior. This must be a column vector.
% Note that a single dummy observation becomes lags+1 rows of the data matrix,
% with a break separating it from the rest of the data. The function treats the
% first lags observations at the top and after each "break" in ydata and xdata as
% initial conditions.
% lambda: weight on "co-persistence" prior dummy observations. This expresses
% belief that when data on *all* y's are stable at their initial levels, they will
% tend to persist at that level. lambda=5 is a reasonable first try. With lambda<0,
% constant term is not included in the dummy observation, so that stationary models
% with means equal to initial ybar do not fit the prior mean. With lambda>0, the prior
% implies that large constants are unlikely if unit roots are present.
% mu: weight on "own persistence" prior dummy observation. Expresses belief
% that when y_i has been stable at its initial level, it will tend to persist
% at that level, regardless of the values of other variables. There is
% one of these for each variable. A reasonable first guess is mu=2.
% The program assumes that the first lags rows of ydata and xdata are real data, not dummies.
% Dummy observations should go at the end, if any. If pre-sample x's are not available,
% repeating the initial xdata(lags+1,:) row or copying xdata(lags+1:2*lags,:) into
% xdata(1:lags,:) are reasonable subsititutes. These values are used in forming the
% persistence priors.
% Code written by Christopher Sims. This version 6/15/03.
[T,nvar]=size(ydata);
nox=isempty(xdata);
if ~nox
[T2,nx]=size(xdata);
else
T2=T;nx=0;xdata=zeros(T2,0);
end
% note that x must be same length as y, even though first part of x will not be used.
% This is so that the lags parameter can be changed without reshaping the xdata matrix.
if T2 ~= T, disp('Mismatch of x and y data lengths'),end
if nargin<4
nbreaks=0;breaks=[];
else
nbreaks=length(breaks);
end
breaks=[0;breaks;T];
smpl=[];
for nb=1:nbreaks+1
smpl=[smpl;[breaks(nb)+lags+1:breaks(nb+1)]'];
end
Tsmpl=size(smpl,1);
X=zeros(Tsmpl,nvar,lags);
for is=1:length(smpl)
X(is,:,:)=ydata(smpl(is)-(1:lags),:)';
end
X=[X(:,:) xdata(smpl,:)];
y=ydata(smpl,:);
% Everything now set up with input data for y=Xb+e
% ------------------Form persistence dummies-------------------
if lambda~=0 | mu>0
ybar=mean(ydata(1:lags,:),1);
if ~nox
xbar=mean(xdata(1:lags,:),1);
else
xbar=[];
end
if lambda~=0
if lambda>0
xdum=lambda*[repmat(ybar,1,lags) xbar];
else
lambda=-lambda;
xdum=lambda*[repmat(ybar,1,lags) zeros(size(xbar))];
end
ydum=zeros(1,nvar);
ydum(1,:)=lambda*ybar;
y=[y;ydum];
X=[X(:,:);xdum];
end
if mu>0
xdum=[repmat(diag(ybar),1,lags) zeros(nvar,nx)]*mu;
ydum=mu*diag(ybar);
X=[X;xdum];
y=[y;ydum];
end
end
[vl,d,vr]=svd(X(:,:),0);
di=1../diag(d);
B=vl'*y;
B=(vr.*repmat(di',nvar*lags+nx,1))*B;
u=y-X(:,:)*B;
xxi=vr.*repmat(di',nvar*lags+nx,1);
xxi=xxi*xxi';
B=reshape(B,[nvar*lags+nx,nvar]); % rhs variables, equations
By=B(1:nvar*lags,:);
By=reshape(By,nvar,lags,nvar);% variables, lags, equations
By=permute(By,[3,1,2]); %equations, variables, lags to match impulsdt.m
if nox
Bx=[];
else
Bx=B(nvar*lags+(1:nx),:)';
end
%logintlh=matrictint(u'*u,xxi,size(X,1)-nvar-1)-.5*nvar*(nvar+1)*log(2*pi);
var.By=By;var.Bx=Bx;var.u=u;var.xxi=xxi;%var.logintlh=logintlh;
% Desired features: 1) automatic dummies for vcv prior
% 2) automatic calculation of integrated pdf, accounting
% for the dummy variables as a prior
% 3) automatic dummies for "Minnesota prior"

View File

@ -1,54 +0,0 @@
function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
%function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
% ydum, xdum: dummy observation data that implement the prior
% breaks: vector of points in the dummy data after which new dummy obs's start
% Set breaks=T+[0;breaks], ydata=[ydata;ydum], xdum=[xdata;xdum], where
% actual data matrix has T rows, in preparing input for rfvar3
% nv,nx,lags: VAR dimensions
% mnprior.tight:Overall tightness of Minnesota prior
% mnprior.decay:Standard deviations of lags shrink as lag^(-decay)
% vprior.sig: Vector of prior modes for diagonal elements of r.f. covariance matrix
% vprior.w: Weight on prior on vcv. 1 corresponds to "one dummy observation" weight
% Should be an integer, and will be rounded if not. vprior.sig is needed
% to scale the Minnesota prior, even if the prior on sigma is not used itself.
% Set vprior.w=0 to achieve this.
% Note: The original Minnesota prior treats own lags asymmetrically, and therefore
% cannot be implemented entirely with dummy observations. It is also usually
% taken to include the sum-of-coefficients and co-persistence components
% that are implemented directly in rfvar3.m. The diagonal prior on v, combined
% with sum-of-coefficients and co-persistence components and with the unit own-first-lag
% prior mean generates larger prior variances for own than for cross-effects even in
% this formulation, but here there is no way to shrink toward a set of unconstrained
% univariate AR's.
%-----------------------
%
if ~isempty(mnprior)
xdum=zeros(lags+1,nx,lags,nv);
ydum=zeros(lags+1,nv,lags,nv);
for il=1:lags
ydum(il+1,:,il,:)=il^mnprior.decay*diag(vprior.sig);
end
ydum(1,:,1,:)=diag(vprior.sig);
ydum=mnprior.tight*reshape(ydum,[lags+1,nv,lags*nv]);
ydum=flipdim(ydum,1);
xdum=mnprior.tight*reshape(xdum,[lags+1,nx,lags*nv]);
xdum=flipdim(xdum,1);
breaks=(lags+1)*[1:(nv*lags)]';
lbreak=breaks(end);
else
ydum=[];
xdum=[];
breaks=[];
lbreak=0;
end
if ~isempty(vprior) & vprior.w>0
ydum2=zeros(lags+1,nv,nv);
xdum2=zeros(lags+1,nx,nv);
ydum2(end,:,:)=diag(vprior.sig);
ydum=cat(3,ydum,ydum2);
xdum=cat(3,xdum,xdum2);
dimy=size(ydum);
ydum=reshape(permute(ydum,[1 3 2]),dimy(1)*dimy(3),nv);
xdum=reshape(permute(xdum,[1 3 2]),dimy(1)*dimy(3),nx);
breaks=[breaks;(lags+1)*[1:nv-1]'+lbreak];
end

View File

@ -815,3 +815,29 @@ PlannerObjectiveStatement::writeOutput(ostream &output, const string &basename)
{
model_tree->writeStaticFile(basename + "_objective");
}
BVARDensityStatement::BVARDensityStatement(int maxnlags_arg, const OptionsList &options_list_arg) :
maxnlags(maxnlags_arg),
options_list(options_list_arg)
{
}
void
BVARDensityStatement::writeOutput(ostream &output, const string &basename) const
{
options_list.writeOutput(output);
output << "bvar_density(" << maxnlags << ");" << endl;
}
BVARForecastStatement::BVARForecastStatement(int nlags_arg, const OptionsList &options_list_arg) :
nlags(nlags_arg),
options_list(options_list_arg)
{
}
void
BVARForecastStatement::writeOutput(ostream &output, const string &basename) const
{
options_list.writeOutput(output);
output << "bvar_forecast(" << nlags << ");" << endl;
}

File diff suppressed because it is too large Load Diff

View File

@ -35,6 +35,10 @@ class ParsingDriver;
%token AR AUTOCORR
%token BAYESIAN_IRF BETA_PDF
%token BVAR_DENSITY BVAR_FORECAST
%token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA
%token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
%token BVAR_REPLIC
%token CALIB CALIB_VAR CHECK CONF_SIG CONSTANT CORR COVAR CUTOFF
%token DATAFILE DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE
%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT
@ -125,6 +129,8 @@ class ParsingDriver;
| model_comparison
| planner_objective
| ramsey_policy
| bvar_density
| bvar_forecast
;
@ -1117,6 +1123,51 @@ markowitz
| o_planner_discount
;
bvar_prior_option : o_bvar_prior_tau
| o_bvar_prior_decay
| o_bvar_prior_lambda
| o_bvar_prior_mu
| o_bvar_prior_omega
| o_bvar_prior_flat
| o_bvar_prior_train
;
bvar_common_option : bvar_prior_option
| o_datafile
| o_xls_sheet
| o_xls_range
| o_first_obs
| o_presample
| o_nobs
| o_prefilter
;
bvar_density_options_list : bvar_common_option COMMA bvar_density_options_list
| bvar_common_option
;
bvar_density : BVAR_DENSITY INT_NUMBER ';'
{ driver.bvar_density($2); }
| BVAR_DENSITY '(' bvar_density_options_list ')' INT_NUMBER ';'
{ driver.bvar_density($5); }
;
bvar_forecast_option : bvar_common_option
| o_forecast
| o_conf_sig
| o_bvar_replic
;
bvar_forecast_options_list : bvar_forecast_option COMMA bvar_forecast_options_list
| bvar_forecast_option
;
bvar_forecast : BVAR_FORECAST INT_NUMBER ';'
{ driver.bvar_forecast($2); }
| BVAR_FORECAST '(' bvar_forecast_options_list ')' INT_NUMBER ';'
{ driver.bvar_forecast($5); }
;
o_dr_algo: DR_ALGO EQUAL INT_NUMBER {driver.option_num("dr_algo", $3);};
o_solve_algo: SOLVE_ALGO EQUAL INT_NUMBER {driver.option_num("solve_algo", $3);};
o_simul_algo: SIMUL_ALGO EQUAL INT_NUMBER {driver.option_num("simul_algo", $3);};
@ -1192,6 +1243,15 @@ markowitz
o_mh_recover : MH_RECOVER {driver.option_num("load_mh_file", "-1");}
o_planner_discount : PLANNER_DISCOUNT EQUAL FLOAT_NUMBER {driver.option_num("planner_discount",$3);}
o_bvar_prior_tau : BVAR_PRIOR_TAU EQUAL signed_float { driver.option_num("bvar_prior_tau", $3); }
o_bvar_prior_decay : BVAR_PRIOR_DECAY EQUAL FLOAT_NUMBER { driver.option_num("bvar_prior_decay", $3); }
o_bvar_prior_lambda : BVAR_PRIOR_LAMBDA EQUAL signed_float { driver.option_num("bvar_prior_lambda", $3); }
o_bvar_prior_mu : BVAR_PRIOR_MU EQUAL FLOAT_NUMBER { driver.option_num("bvar_prior_mu", $3); }
o_bvar_prior_omega : BVAR_PRIOR_OMEGA EQUAL INT_NUMBER { driver.option_num("bvar_prior_omega", $3); }
o_bvar_prior_flat : BVAR_PRIOR_FLAT { driver.option_num("bvar_prior_flat", "1"); }
o_bvar_prior_train : BVAR_PRIOR_TRAIN EQUAL INT_NUMBER { driver.option_num("bvar_prior_train", $3); }
o_bvar_replic : BVAR_REPLIC EQUAL INT_NUMBER { driver.option_num("bvar_replic", $3); }
range : NAME ':' NAME
{
$1->append(":");

View File

@ -83,6 +83,9 @@ int sigma_e = 0;
<INITIAL>planner_objective {BEGIN DYNARE_STATEMENT; return token::PLANNER_OBJECTIVE;}
<INITIAL>ramsey_policy {BEGIN DYNARE_STATEMENT; return token::RAMSEY_POLICY;}
<INITIAL>bvar_density {BEGIN DYNARE_STATEMENT; return token::BVAR_DENSITY; }
<INITIAL>bvar_forecast {BEGIN DYNARE_STATEMENT; return token::BVAR_FORECAST; }
/* End of a Dynare statement */
<DYNARE_STATEMENT>; {
if (!sigma_e)
@ -159,6 +162,15 @@ int sigma_e = 0;
<DYNARE_STATEMENT>noconstant {return token::NOCONSTANT;}
<DYNARE_STATEMENT>covar {return token::COVAR;}
<DYNARE_STATEMENT>bvar_prior_tau { return token::BVAR_PRIOR_TAU; }
<DYNARE_STATEMENT>bvar_prior_decay { return token::BVAR_PRIOR_DECAY; }
<DYNARE_STATEMENT>bvar_prior_lambda { return token::BVAR_PRIOR_LAMBDA; }
<DYNARE_STATEMENT>bvar_prior_mu { return token::BVAR_PRIOR_MU; }
<DYNARE_STATEMENT>bvar_prior_omega { return token::BVAR_PRIOR_OMEGA; }
<DYNARE_STATEMENT>bvar_prior_flat { return token::BVAR_PRIOR_FLAT; }
<DYNARE_STATEMENT>bvar_prior_train { return token::BVAR_PRIOR_TRAIN; }
<DYNARE_STATEMENT>bvar_replic { return token::BVAR_REPLIC; }
<DYNARE_STATEMENT>[\$][^$]*[\$] {
strtok(yytext+1, "$");
yylval->string_val = new string(yytext + 1);

View File

@ -992,6 +992,22 @@ ParsingDriver::ramsey_policy()
options_list.clear();
}
void
ParsingDriver::bvar_density(string *maxnlags)
{
mod_file->addStatement(new BVARDensityStatement(atoi(maxnlags->c_str()), options_list));
options_list.clear();
delete maxnlags;
}
void
ParsingDriver::bvar_forecast(string *nlags)
{
mod_file->addStatement(new BVARForecastStatement(atoi(nlags->c_str()), options_list));
options_list.clear();
delete nlags;
}
NodeID
ParsingDriver::add_model_equal(NodeID arg1, NodeID arg2)
{

View File

@ -403,4 +403,25 @@ public:
virtual void writeOutput(ostream &output, const string &basename) const;
};
class BVARDensityStatement : public Statement
{
private:
const int maxnlags;
const OptionsList options_list;
public:
BVARDensityStatement(int maxnlags_arg, const OptionsList &options_list_arg);
virtual void writeOutput(ostream &output, const string &basename) const;
};
class BVARForecastStatement : public Statement
{
private:
const int nlags;
const OptionsList options_list;
public:
BVARForecastStatement(int nlags_arg, const OptionsList &options_list_arg);
virtual void writeOutput(ostream &output, const string &basename) const;
};
#endif

View File

@ -131,149 +131,159 @@ namespace yy
AUTOCORR = 259,
BAYESIAN_IRF = 260,
BETA_PDF = 261,
CALIB = 262,
CALIB_VAR = 263,
CHECK = 264,
CONF_SIG = 265,
CONSTANT = 266,
CORR = 267,
COVAR = 268,
CUTOFF = 269,
DATAFILE = 270,
DR_ALGO = 271,
DROP = 272,
DSAMPLE = 273,
DYNASAVE = 274,
DYNATYPE = 275,
END = 276,
ENDVAL = 277,
EQUAL = 278,
ESTIMATION = 279,
ESTIMATED_PARAMS = 280,
ESTIMATED_PARAMS_BOUNDS = 281,
ESTIMATED_PARAMS_INIT = 282,
FILENAME = 283,
FILTER_STEP_AHEAD = 284,
FILTERED_VARS = 285,
FIRST_OBS = 286,
FLOAT_NUMBER = 287,
FORECAST = 288,
GAMMA_PDF = 289,
GCC_COMPILER = 290,
GRAPH = 291,
HISTVAL = 292,
HP_FILTER = 293,
HP_NGRID = 294,
INITVAL = 295,
INT_NUMBER = 296,
INV_GAMMA_PDF = 297,
IRF = 298,
KALMAN_ALGO = 299,
KALMAN_TOL = 300,
LAPLACE = 301,
LCC_COMPILER = 302,
LIK_ALGO = 303,
LIK_INIT = 304,
LINEAR = 305,
LOAD_MH_FILE = 306,
LOGLINEAR = 307,
MARKOWITZ = 308,
MH_DROP = 309,
MH_INIT_SCALE = 310,
MH_JSCALE = 311,
MH_MODE = 312,
MH_NBLOCKS = 313,
MH_REPLIC = 314,
MH_RECOVER = 315,
MODE_CHECK = 316,
MODE_COMPUTE = 317,
MODE_FILE = 318,
MODEL = 319,
MODEL_COMPARISON = 320,
MSHOCKS = 321,
MODEL_COMPARISON_APPROXIMATION = 322,
MODIFIEDHARMONICMEAN = 323,
MOMENTS_VARENDO = 324,
NAME = 325,
NOBS = 326,
NOCONSTANT = 327,
NOCORR = 328,
NODIAGNOSTIC = 329,
NOFUNCTIONS = 330,
NOGRAPH = 331,
NOMOMENTS = 332,
NOPRINT = 333,
NORMAL_PDF = 334,
OBSERVATION_TRENDS = 335,
OLR = 336,
OLR_INST = 337,
OLR_BETA = 338,
OPTIM = 339,
OPTIM_WEIGHTS = 340,
ORDER = 341,
OSR = 342,
OSR_PARAMS = 343,
PARAMETERS = 344,
PERIODS = 345,
PLANNER_OBJECTIVE = 346,
PREFILTER = 347,
PRESAMPLE = 348,
PRINT = 349,
PRIOR_TRUNC = 350,
PRIOR_ANALYSIS = 351,
POSTERIOR_ANALYSIS = 352,
QZ_CRITERIUM = 353,
RELATIVE_IRF = 354,
REPLIC = 355,
RPLOT = 356,
SHOCKS = 357,
SIGMA_E = 358,
SIMUL = 359,
SIMUL_ALGO = 360,
SIMUL_SEED = 361,
SMOOTHER = 362,
SOLVE_ALGO = 363,
SPARSE_DLL = 364,
STDERR = 365,
STEADY = 366,
STOCH_SIMUL = 367,
TEX = 368,
RAMSEY_POLICY = 369,
PLANNER_DISCOUNT = 370,
TEX_NAME = 371,
UNIFORM_PDF = 372,
UNIT_ROOT_VARS = 373,
USE_DLL = 374,
VALUES = 375,
VAR = 376,
VAREXO = 377,
VAREXO_DET = 378,
VAROBS = 379,
XLS_SHEET = 380,
XLS_RANGE = 381,
COMMA = 382,
MINUS = 383,
PLUS = 384,
DIVIDE = 385,
TIMES = 386,
UMINUS = 387,
POWER = 388,
EXP = 389,
LOG = 390,
LOG10 = 391,
SIN = 392,
COS = 393,
TAN = 394,
ASIN = 395,
ACOS = 396,
ATAN = 397,
SINH = 398,
COSH = 399,
TANH = 400,
ASINH = 401,
ACOSH = 402,
ATANH = 403,
SQRT = 404
BVAR_DENSITY = 262,
BVAR_FORECAST = 263,
BVAR_PRIOR_DECAY = 264,
BVAR_PRIOR_FLAT = 265,
BVAR_PRIOR_LAMBDA = 266,
BVAR_PRIOR_MU = 267,
BVAR_PRIOR_OMEGA = 268,
BVAR_PRIOR_TAU = 269,
BVAR_PRIOR_TRAIN = 270,
BVAR_REPLIC = 271,
CALIB = 272,
CALIB_VAR = 273,
CHECK = 274,
CONF_SIG = 275,
CONSTANT = 276,
CORR = 277,
COVAR = 278,
CUTOFF = 279,
DATAFILE = 280,
DR_ALGO = 281,
DROP = 282,
DSAMPLE = 283,
DYNASAVE = 284,
DYNATYPE = 285,
END = 286,
ENDVAL = 287,
EQUAL = 288,
ESTIMATION = 289,
ESTIMATED_PARAMS = 290,
ESTIMATED_PARAMS_BOUNDS = 291,
ESTIMATED_PARAMS_INIT = 292,
FILENAME = 293,
FILTER_STEP_AHEAD = 294,
FILTERED_VARS = 295,
FIRST_OBS = 296,
FLOAT_NUMBER = 297,
FORECAST = 298,
GAMMA_PDF = 299,
GCC_COMPILER = 300,
GRAPH = 301,
HISTVAL = 302,
HP_FILTER = 303,
HP_NGRID = 304,
INITVAL = 305,
INT_NUMBER = 306,
INV_GAMMA_PDF = 307,
IRF = 308,
KALMAN_ALGO = 309,
KALMAN_TOL = 310,
LAPLACE = 311,
LCC_COMPILER = 312,
LIK_ALGO = 313,
LIK_INIT = 314,
LINEAR = 315,
LOAD_MH_FILE = 316,
LOGLINEAR = 317,
MARKOWITZ = 318,
MH_DROP = 319,
MH_INIT_SCALE = 320,
MH_JSCALE = 321,
MH_MODE = 322,
MH_NBLOCKS = 323,
MH_REPLIC = 324,
MH_RECOVER = 325,
MODE_CHECK = 326,
MODE_COMPUTE = 327,
MODE_FILE = 328,
MODEL = 329,
MODEL_COMPARISON = 330,
MSHOCKS = 331,
MODEL_COMPARISON_APPROXIMATION = 332,
MODIFIEDHARMONICMEAN = 333,
MOMENTS_VARENDO = 334,
NAME = 335,
NOBS = 336,
NOCONSTANT = 337,
NOCORR = 338,
NODIAGNOSTIC = 339,
NOFUNCTIONS = 340,
NOGRAPH = 341,
NOMOMENTS = 342,
NOPRINT = 343,
NORMAL_PDF = 344,
OBSERVATION_TRENDS = 345,
OLR = 346,
OLR_INST = 347,
OLR_BETA = 348,
OPTIM = 349,
OPTIM_WEIGHTS = 350,
ORDER = 351,
OSR = 352,
OSR_PARAMS = 353,
PARAMETERS = 354,
PERIODS = 355,
PLANNER_OBJECTIVE = 356,
PREFILTER = 357,
PRESAMPLE = 358,
PRINT = 359,
PRIOR_TRUNC = 360,
PRIOR_ANALYSIS = 361,
POSTERIOR_ANALYSIS = 362,
QZ_CRITERIUM = 363,
RELATIVE_IRF = 364,
REPLIC = 365,
RPLOT = 366,
SHOCKS = 367,
SIGMA_E = 368,
SIMUL = 369,
SIMUL_ALGO = 370,
SIMUL_SEED = 371,
SMOOTHER = 372,
SOLVE_ALGO = 373,
SPARSE_DLL = 374,
STDERR = 375,
STEADY = 376,
STOCH_SIMUL = 377,
TEX = 378,
RAMSEY_POLICY = 379,
PLANNER_DISCOUNT = 380,
TEX_NAME = 381,
UNIFORM_PDF = 382,
UNIT_ROOT_VARS = 383,
USE_DLL = 384,
VALUES = 385,
VAR = 386,
VAREXO = 387,
VAREXO_DET = 388,
VAROBS = 389,
XLS_SHEET = 390,
XLS_RANGE = 391,
COMMA = 392,
MINUS = 393,
PLUS = 394,
DIVIDE = 395,
TIMES = 396,
UMINUS = 397,
POWER = 398,
EXP = 399,
LOG = 400,
LOG10 = 401,
SIN = 402,
COS = 403,
TAN = 404,
ASIN = 405,
ACOS = 406,
ATAN = 407,
SINH = 408,
COSH = 409,
TANH = 410,
ASINH = 411,
ACOSH = 412,
ATANH = 413,
SQRT = 414
};
};

View File

@ -306,6 +306,10 @@ public:
void end_planner_objective(NodeID expr);
//! ramsey policy statement
void ramsey_policy();
//! BVAR marginal density
void bvar_density(string *maxnlags);
//! BVAR forecast
void bvar_forecast(string *nlags);
//! Writes token "arg1=arg2" to model tree
NodeID add_model_equal(NodeID arg1, NodeID arg2);
//! Writes token "arg=0" to model tree