From dc7fca7ece954fab0932d075305b2d201aae74d3 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 24 Nov 2016 12:34:52 +0100 Subject: [PATCH] var_forecast: use example1 in forecast, add code to use estimation via rfvar3 --- matlab/bvar_toolbox.m | 105 ---------------------------- matlab/rfvar3.m | 122 +++++++++++++++++++++++++++++++++ tests/ECB/example1.mod | 51 ++++++++++++++ tests/ECB/example1_var.mod | 53 ++++++++++++++ tests/ECB/run_var_comparison.m | 49 +++++++++---- 5 files changed, 262 insertions(+), 118 deletions(-) create mode 100644 matlab/rfvar3.m create mode 100644 tests/ECB/example1.mod create mode 100644 tests/ECB/example1_var.mod diff --git a/matlab/bvar_toolbox.m b/matlab/bvar_toolbox.m index 76c95ef9d..05086a875 100644 --- a/matlab/bvar_toolbox.m +++ b/matlab/bvar_toolbox.m @@ -216,108 +216,3 @@ 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(1:(end-1)); - - -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. - -% Original file downloaded from: -% http://sims.princeton.edu/yftp/VARtools/matlab/rfvar3.m - -[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; diff --git a/matlab/rfvar3.m b/matlab/rfvar3.m new file mode 100644 index 000000000..25c5ce8d1 --- /dev/null +++ b/matlab/rfvar3.m @@ -0,0 +1,122 @@ +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. + +% Original file downloaded from: +% http://sims.princeton.edu/yftp/VARtools/matlab/rfvar3.m + +% Copyright (C) 2003-2007 Christopher Sims +% Copyright (C) 2007-2012 Dynare Team +% +% 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 . + +[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; +end \ No newline at end of file diff --git a/tests/ECB/example1.mod b/tests/ECB/example1.mod new file mode 100644 index 000000000..8aa4b7343 --- /dev/null +++ b/tests/ECB/example1.mod @@ -0,0 +1,51 @@ +// Example 1 from Collard's guide to Dynare +var y, c, k, a, h, b; +varexo e, u; + +verbatim; +% I want these comments included in +% example1.m 1999q1 1999y +% +var = 1; +end; + +parameters beta, rho, alpha, delta, theta, psi, tau; + +alpha = 0.36; +rho = 0.95; +tau = 0.025; +beta = 0.99; +delta = 0.025; +psi = 0; +theta = 2.95; + +phi = 0.1; + +model; +c*theta*h^(1+psi)=(1-alpha)*y; +k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1))) + *(exp(b(+1))*alpha*y(+1)+(1-delta)*k)); +y = exp(a)*(k(-1)^alpha)*(h^(1-alpha)); +k = exp(b)*(y-c)+(1-delta)*k(-1); +a = rho*a(-1)+tau*b(-1) + e; +b = tau*a(-1)+rho*b(-1) + u; +end; + +initval; +y = 1.08068253095672; +c = 0.80359242014163; +h = 0.29175631001732; +k = 11.08360443260358; +a = 0; +b = 0; +e = 0; +u = 0; +end; + +shocks; +var e; stderr 0.009; +var u; stderr 0.009; +var e, u = phi*0.009*0.009; +end; + +stoch_simul(order=1, periods=200); diff --git a/tests/ECB/example1_var.mod b/tests/ECB/example1_var.mod new file mode 100644 index 000000000..b1d705548 --- /dev/null +++ b/tests/ECB/example1_var.mod @@ -0,0 +1,53 @@ +// Example 1 from Collard's guide to Dynare +var y, c, k, a, h, b; +varexo e, u; + +verbatim; +% I want these comments included in +% example1.m 1999q1 1999y +% +var = 1; +end; + +parameters beta, rho, alpha, delta, theta, psi, tau; + +alpha = 0.36; +rho = 0.95; +tau = 0.025; +beta = 0.99; +delta = 0.025; +psi = 0; +theta = 2.95; + +phi = 0.1; + +var_model(model_name=my_var_est, order=1) y c; + +model; +c*theta*h^(1+psi)=(1-alpha)*y; +k = beta*(((exp(b)*c)/(exp(b(+1))*c(1))) + *(exp(b(+1))*alpha*var_expectation(y, model_name=my_var_est)+(1-delta)*k)); +y = exp(a)*(k(-1)^alpha)*(h^(1-alpha)); +k = exp(b)*(y-c)+(1-delta)*k(-1); +a = rho*a(-1)+tau*b(-1) + e; +b = tau*a(-1)+rho*b(-1) + u; +end; + +initval; +y = 1.08068253095672; +c = 0.80359242014163; +h = 0.29175631001732; +k = 11.08360443260358; +a = 0; +b = 0; +e = 0; +u = 0; +end; + +shocks; +var e; stderr 0.009; +var u; stderr 0.009; +var e, u = phi*0.009*0.009; +end; + +stoch_simul(order=1, periods=200); diff --git a/tests/ECB/run_var_comparison.m b/tests/ECB/run_var_comparison.m index 5677b5d09..942ffe92e 100644 --- a/tests/ECB/run_var_comparison.m +++ b/tests/ECB/run_var_comparison.m @@ -1,30 +1,45 @@ clearvars clearvars -global +close all + !rm nkm_saved_data.mat !rm my_var_est.mat !rm nkm_var_saved_data.mat %% Call Dynare without VAR -dynare nkm.mod +dynare example1.mod save('nkm_saved_data.mat') %% Estimate VAR using simulated data disp('VAR Estimation'); -Y = oo_.endo_simul(2:3, 2:end); -Z = [ones(1, size(Y,2)); oo_.endo_simul(2:3, 1:end-1)]; -% OLS -B = Y*transpose(Z)/(Z*transpose(Z)); -%% save values +% MLS estimate of mu and B (autoregressive_matrices) +% Y = mu + B*Z +% from New Introduction to Multiple Time Series Analysis +Y = oo_.endo_simul(1:2, 2:end); +Z = [ ... + ones(1, size(Y,2)); ... + oo_.endo_simul(1:2, 1:end-1); ... + ]; +%B = Y*Z'*inv(Z*Z'); +B = Y*Z'/(Z*Z'); mu = B(:, 1); autoregressive_matrices{1} = B(:, 2:end); + +% Sims +% (provides same result as above) +% var = rfvar3(Y',1,zeros(size(Y')),0,5,2) +% mu = var.B(:, 1); +% autoregressive_matrices{1} = var.B(:, 2:end); + +%% save values save('my_var_est.mat', 'mu', 'autoregressive_matrices'); %% Call Dynare with VAR clearvars clearvars -global -dynare nkm_var.mod +dynare example1_var.mod save('nkm_var_saved_data.mat') %% compare values @@ -35,11 +50,19 @@ zerotol = 1e-12; nv = load('nkm_saved_data.mat'); wv = load('nkm_var_saved_data.mat'); -assert(max(max(abs(nv.y - wv.y))) < zerotol); -assert(max(max(abs(nv.pi - wv.pi))) < zerotol); -assert(max(max(abs(nv.i - wv.i))) < zerotol); +ridx = 3; +cidx = 2; +exo_names = nv.M_.exo_names; +endo_names = nv.M_.endo_names; -fn = fieldnames(nv.oo_.irfs); -for i=1:length(fn) - assert(max(max(abs(nv.oo_.irfs.(fn{i}) - (wv.oo_.irfs.(fn{i}))))) < zerotol); +for i = 1:length(exo_names) + figure('Name', ['Shock to ' exo_names(i)]); + for j = 1:length(endo_names) + subplot(ridx, cidx, j); + hold on + title(endo_names(j)); + plot(nv.oo_.irfs.([endo_names(j) '_' exo_names(i)])); + plot(wv.oo_.irfs.([endo_names(j) '_' exo_names(i)]), '--'); + hold off + end end \ No newline at end of file