From 683c89de43b44811e9135abe62cf5c7db7facb48 Mon Sep 17 00:00:00 2001
From: Houtan Bastani <houtan@dynare.org>
Date: Wed, 24 Jan 2018 15:22:30 +0100
Subject: [PATCH] SUR: add option to estimate certain parameters

---
 matlab/ols/sur.m  | 41 ++++++++++++++++++++++++++++++++++-------
 matlab/surgibbs.m | 22 ++++------------------
 2 files changed, 38 insertions(+), 25 deletions(-)

diff --git a/matlab/ols/sur.m b/matlab/ols/sur.m
index b68df21e8..eba7069e8 100644
--- a/matlab/ols/sur.m
+++ b/matlab/ols/sur.m
@@ -1,9 +1,10 @@
-function varargout = sur(ds, eqtags)
-% function varargout = sur(ds, eqtags)
+function varargout = sur(ds, param_names, eqtags)
+%function varargout = sur(ds, param_names, eqtags)
 % Seemingly Unrelated Regressions
 %
 % INPUTS
 %   ds                [dseries]    data to use in estimation
+%   param_names       [cellstr]    list of parameters to estimate
 %   eqtags            [cellstr]    names of equation tags to estimate. If empty,
 %                                  estimate all equations
 %
@@ -33,8 +34,13 @@ function varargout = sur(ds, eqtags)
 global M_ oo_ options_
 
 %% Check input argument
-assert(nargin == 1 || nargin == 2, 'You must provide one or two arguments');
+assert(nargin >= 1 && nargin <= 3, 'You must provide one, two, or three arguments');
 assert(~isempty(ds) && isdseries(ds), 'The first argument must be a dseries');
+if nargin >= 2
+    assert(iscellstr(param_names), 'The 2nd argument must be a cellstr');
+else
+    param_names = {};
+end
 
 %% Read JSON
 jsonfile = [M_.fname '_original.json'];
@@ -44,12 +50,12 @@ end
 
 jsonmodel = loadjson(jsonfile);
 jsonmodel = jsonmodel.model;
-if nargin == 2
+if nargin == 3
     jsonmodel = getEquationsByTags(jsonmodel, 'name', eqtags);
 end
 
 %% Find parameters and variable names in equations and setup estimation matrices
-[X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, pidxs, surconstrainedparams] = ...
+[X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, opidxs, surconstrainedparams] = ...
     pooled_sur_common(ds, jsonmodel);
 
 if nargin == 1 && size(X, 2) ~= M_.param_nbr
@@ -79,11 +85,32 @@ for i = 1:length(jsonmodel)
     end
 end
 
+if ~isempty(param_names)
+    pnamesall = M_.param_names(opidxs);
+    nparams = length(param_names);
+    pidxs = zeros(nparams, 1);
+    for i = 1:nparams
+        idxs = find(strcmp(param_names{i}, pnamesall));
+        if isempty(idxs)
+            if ~isempty(eqtags)
+                error(['Could not find ' param_names{i} ...
+                    ' in the provided equations specified by ' strjoin(eqtags, ',')]);
+            end
+            error('Unspecified error. Please report');
+        end
+        pidxs(i) = idxs;
+    end
+    vars = [vars{:}];
+    vars = {vars(pidxs)};
+    newY = newY - newX(:, setdiff(1:size(newX, 2), pidxs)) * M_.params(setdiff(opidxs, opidxs(pidxs), 'stable'));
+    newX = newX(:, pidxs);
+end
+
 %% Return to surgibbs if called from there
 st = dbstack(1);
 if strcmp(st(1).name, 'surgibbs')
     varargout{1} = length(maxfp:minlp); %dof
-    varargout{2} = pidxs;
+    varargout{2} = opidxs(pidxs);
     varargout{3} = newX;
     varargout{4} = newY;
     varargout{5} = length(jsonmodel);
@@ -106,7 +133,7 @@ kLeye = kron(chol(inv(M_.Sigma_e)), eye(oo_.sur.dof));
 [q, r] = qr(kLeye*X, 0);
 oo_.sur.beta = r\(q'*kLeye*Y);
 
-M_.params(pidxs, 1) = oo_.sur.beta;
+M_.params(opidxs(pidxs)) = oo_.sur.beta;
 
 % Yhat
 oo_.sur.Yhat = X * oo_.sur.beta;
diff --git a/matlab/surgibbs.m b/matlab/surgibbs.m
index 13ed43501..3df6e1533 100644
--- a/matlab/surgibbs.m
+++ b/matlab/surgibbs.m
@@ -70,29 +70,15 @@ end
 
 %% Estimation
 if nargin == 8
-    [nobs, pidxs, X, Y, m] = sur(ds, eqtags);
+    [nobs, pidxs, X, Y, m] = sur(ds, param_names, eqtags);
 else
-    [nobs, pidxs, X, Y, m] = sur(ds);
+    [nobs, pidxs, X, Y, m] = sur(ds, param_names);
 end
-pnamesall = M_.param_names(pidxs);
-nparams = length(param_names);
-pidxs = zeros(nparams, 1);
-for i = 1:nparams
-    idxs = find(strcmp(param_names{i}, pnamesall));
-    if isempty(idxs)
-        if ~isempty(eqtags)
-            error(['Could not find ' param_names{i} ...
-                ' in the provided equations specified by ' strjoin(eqtags, ',')]);
-        end
-        error('Unspecified error. Please report');
-    end
-    pidxs(i) = idxs;
-end
-X = X(:, pidxs);
 beta = beta0;
 A = inv(A);
 thinidx = 1;
 drawidx = 1;
+nparams = length(param_names);
 oo_.surgibbs.betadraws = zeros(floor((ndraws-discarddraws)/thin), nparams);
 for i = 1:ndraws
     % Draw Omega, given X, Y, Beta
@@ -119,7 +105,7 @@ end
 
 % save parameter values
 oo_.surgibbs.beta = (sum(oo_.surgibbs.betadraws)/rows(oo_.surgibbs.betadraws))';
-M_.params(pidxs, 1) = oo_.surgibbs.beta;
+M_.params(pidxs) = oo_.surgibbs.beta;
 
 %% Print Output
 dyn_table('Gibbs Sampling on SUR', {}, {}, param_names, ...