From 12cd0d819438ce8aa8a91bcd12a23b9903a1f38f Mon Sep 17 00:00:00 2001
From: Houtan Bastani <houtan@dynare.org>
Date: Fri, 7 Jul 2017 14:22:52 +0200
Subject: [PATCH] ols: add test statistics

---
 matlab/olseqs.m                      | 85 +++++++++++++++++++++++++---
 tests/ecb/contribution-plots/var.mod |  2 +
 2 files changed, 80 insertions(+), 7 deletions(-)

diff --git a/matlab/olseqs.m b/matlab/olseqs.m
index 283332fad..a263a08af 100644
--- a/matlab/olseqs.m
+++ b/matlab/olseqs.m
@@ -1,12 +1,12 @@
-function B = olseqs(ds, varargin)
-%function B = olseqs(ds, varargin)
+function olseqs(ds, varargin)
+%function olseqs(ds, varargin)
 % Run OLS on chosen model equations
 %
 % INPUTS
 %   ds      [dseries]    data
 %
 % OUTPUTS
-%   B       [vector]     estimated parameters
+%   none
 %
 % SPECIAL REQUIREMENTS
 %   none
@@ -28,8 +28,9 @@ function B = olseqs(ds, varargin)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global M_
+global M_ oo_
 
+%% Check input
 assert(nargin == 1 || nargin == 3, 'Incorrect number of arguments passed to olseqs');
 
 jsonfile = [M_.fname '_original.json'];
@@ -37,15 +38,18 @@ if exist(jsonfile, 'file') ~= 2
     error('Could not find %s! Please use the json option (See the Dynare invocation section in the reference manual).', jsonfile);
 end
 
+%% Get Equation(s)
 jsonmodel = loadjson(jsonfile);
 jsonmodel = jsonmodel.model;
 [lhs, rhs, lineno] = getEquationsByTags(jsonmodel, varargin{:});
 
+%% Construct regression matrices
 Y = ds{lhs}.data;
 
 rhs_ = strsplit(rhs, {'+','-','*','/','^','log(','exp(','(',')'});
 rhs_(cellfun(@(x) all(isstrprop(x, 'digit')), rhs_)) = [];
 vnames = setdiff(rhs_, cellstr(M_.param_names));
+pnames = setdiff(rhs_, vnames);
 regexprnoleads = cell2mat(strcat('(', vnames, {'\(\d+\))|'}));
 if ~isempty(regexp(rhs, regexprnoleads(1:end-1), 'match'))
     error(['olseqs: you cannot have leads in equation on line ' lineno ': ' lhs ' = ' rhs]);
@@ -59,7 +63,74 @@ X = cell2mat(cellfun(@eval, strcat('ds.', vwlags, '.data'), 'UniformOutput', fal
 Y = Y(row+1:end, :);
 X = X(row+1:end, :);
 
-% OLS Estimation
-B = (X'*X)\X'*Y;
+% Add intercept
+% X = [ones(size(X,1), 1), X];
 
-end
\ No newline at end of file
+%% OLS Estimation
+% From LeSage, James P. "Applied Econometrics using MATLAB"
+tagv = varargin{2};
+[nobs, nvars] = size(X);
+oo_.ols.(tagv).dof = nobs - nvars;
+
+% Estimated Parameters
+[q, r] = qr(X, 0);
+xpxi = (r'*r)\eye(nvars);
+oo_.ols.(tagv).beta = r\(q'*Y);
+
+% Yhat
+oo_.ols.(tagv).Yhat = X*oo_.ols.(tagv).beta;
+
+% Residuals
+oo_.ols.(tagv).resid = Y - oo_.ols.(tagv).Yhat;
+
+% Estimate for sigma^2
+SS_res = oo_.ols.(tagv).resid'*oo_.ols.(tagv).resid;
+oo_.ols.(tagv).s2 = SS_res/oo_.ols.(tagv).dof;
+
+% R^2
+ym = Y - mean(Y);
+SS_tot = ym'*ym;
+oo_.ols.(tagv).R2 = 1 - SS_res/SS_tot;
+
+% Adjusted R^2
+oo_.ols.(tagv).adjR2 = oo_.ols.(tagv).R2 - (1 - oo_.ols.(tagv).R2)*nvars/(oo_.ols.(tagv).dof-1);
+
+% Durbin-Watson
+ediff = oo_.ols.(tagv).resid(2:nobs) - oo_.ols.(tagv).resid(1:nobs-1);
+oo_.ols.(tagv).dw = (ediff'*ediff)/SS_res;
+
+% Standard Error
+oo_.ols.(tagv).stderr = sqrt(oo_.ols.(tagv).s2*diag(xpxi));
+
+% T-Stat
+oo_.ols.(tagv).tstat = oo_.ols.(tagv).beta./oo_.ols.(tagv).stderr;
+
+%% Print Output
+fprintf('OLS Estimation of equation on line %d of %s\n', lineno, [M_.fname '.mod']);
+fprintf('Dependent Variable: %s\n', lhs);
+fprintf('No. Independent Variables: %d\n', nvars);
+fprintf('Observations: %d\n', nobs);
+maxstrlen = 0;
+for i=1:length(vwlags)
+    slen = length(vwlags{i});
+    if slen > maxstrlen
+        maxstrlen = slen;
+    end
+end
+titlespacing = repmat(' ', 1, 4 + maxstrlen + 4) ;
+fprintf('%sCoefficients    t-statistic      Std. Error\n', titlespacing);
+fprintf('%s____________    ____________    ____________\n\n', titlespacing);
+format = ['    %-' num2str(maxstrlen) 's'];
+for i = 1:length(vwlags)
+    fprintf(format, vwlags{i});
+    fprintf('%12.5f    %12.5f    %12.5f\n', ...
+        oo_.ols.(tagv).beta(i), ...
+        oo_.ols.(tagv).tstat(i), ...
+        oo_.ols.(tagv).stderr(i));
+end
+fprintf('\nR^2: %f\n', oo_.ols.(tagv).R2);
+fprintf('R^2 Adjusted: %f\n', oo_.ols.(tagv).adjR2);
+fprintf('s^2: %f\n', oo_.ols.(tagv).s2);
+fprintf('Durbin-Watson: %f\n', oo_.ols.(tagv).dw);
+fprintf('%s\n', repmat('-', 1, 4 + maxstrlen + 4 + 44));
+end
diff --git a/tests/ecb/contribution-plots/var.mod b/tests/ecb/contribution-plots/var.mod
index dc8c300a6..18fec8119 100644
--- a/tests/ecb/contribution-plots/var.mod
+++ b/tests/ecb/contribution-plots/var.mod
@@ -28,4 +28,6 @@ ds1 = dseries(randn(30, 3), 1, {'ffr', 'unrate', 'cpi'});
 // Baseline paths for the variables.
 ds0 = dseries(zeros(30, 3), 1, {'ffr', 'unrate', 'cpi'});
 
+olseqs(ds1, 'eqnum', 'ffr');
+
 plot_contributions('eqnum', 'ffr', ds1, ds0);
\ No newline at end of file