From cd7658541a9615e050d3afd23ab1d978e866fecd Mon Sep 17 00:00:00 2001
From: Houtan <houtanb@gmail.com>
Date: Tue, 12 Dec 2017 13:07:25 +0100
Subject: [PATCH] sur: handle case of constrained parameter

---
 matlab/ols/pooled_ols.m        |  2 +-
 matlab/ols/pooled_sur_common.m | 14 +++++++++++---
 matlab/ols/sur.m               | 10 ++++++++--
 3 files changed, 20 insertions(+), 6 deletions(-)

diff --git a/matlab/ols/pooled_ols.m b/matlab/ols/pooled_ols.m
index 38ef18255..1ec12b70a 100644
--- a/matlab/ols/pooled_ols.m
+++ b/matlab/ols/pooled_ols.m
@@ -87,7 +87,7 @@ end
 %% Find parameters and variable names in every equation & Setup estimation matrices
 M_exo_names_trim = cellstr(M_.exo_names);
 M_param_names_trim = cellstr(M_.param_names);
-[X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, pidxs] = ...
+[X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, surpidxs, surconstrainedparams] = ...
     pooled_sur_common(ds, lhs, rhs, lineno, M_exo_names_trim, M_param_names_trim);
 
 if overlapping_dates
diff --git a/matlab/ols/pooled_sur_common.m b/matlab/ols/pooled_sur_common.m
index a63272e83..9548c5847 100644
--- a/matlab/ols/pooled_sur_common.m
+++ b/matlab/ols/pooled_sur_common.m
@@ -1,4 +1,4 @@
-function [X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, surpidxs] = pooled_sur_common(ds, lhs, rhs, lineno, M_exo_names_trim, M_param_names_trim);
+function [X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, surpidxs, surconstrainedparams] = pooled_sur_common(ds, lhs, rhs, lineno, M_exo_names_trim, M_param_names_trim);
 
 %
 % Code common to sur.m and pooled_ols.m
@@ -45,6 +45,7 @@ enddates = cell(length(lhs), 1);
 residnames = cell(length(lhs), 1);
 surpidxs = zeros(M_.param_nbr, 1); 
 surpidx = 0; 
+surconstrainedparams = [];
 for i = 1:length(lhs)
     rhs_ = strsplit(rhs{i}, {'+','-','*','/','^','log(','ln(','log10(','exp(','(',')','diff('});
     rhs_(cellfun(@(x) all(isstrprop(x, 'digit')), rhs_)) = [];
@@ -61,17 +62,20 @@ for i = 1:length(lhs)
     vnames = cell(1, length(pnames));
     splitstrings = cell(length(pnames), 1);
     xjdata = dseries;
+    dropvname = true(1,length(pnames));
     for j = 1:length(pnames)
         createdvar = false;
         idx = find(strcmp(pbeta, pnames{j}));
         if isempty(idx)
             pbeta = [pbeta; pnames{j}];
             pidxs(j) = length(pbeta);
+            surpidx = surpidx + 1;
+            surpidxs(surpidx, 1) = find(strcmp(pnames{j}, M_param_names_trim));
         else
             pidxs(j) = idx;
+            surconstrainedparams = [surconstrainedparams idx];
+            dropvname(j) = false;
         end
-        surpidx = surpidx + 1;
-        surpidxs(surpidx, 1) = find(strcmp(pnames{j}, M_param_names_trim));
 
         pregex = [...
             mathops pnames{j} mathops ...
@@ -119,6 +123,9 @@ for i = 1:length(lhs)
         xjdatatmp.rename_(num2str(j));
         xjdata = [xjdata xjdatatmp];
     end
+    if ~all(dropvname)
+        vnames = vnames(dropvname);
+    end
 
     lhssub = getRhsToSubFromLhs(ds, rhs{i}, regex, [splitstrings; pnames]);
 
@@ -154,4 +161,5 @@ for i = 1:length(lhs)
         X(startidxs(i):startidxs(i)+lp-fp, pidxs) = xjdata(fp:lp).data;
     end
 end
+surpidxs = surpidxs(1:surpidx);
 end
\ No newline at end of file
diff --git a/matlab/ols/sur.m b/matlab/ols/sur.m
index 75becfe1a..4913c35c7 100644
--- a/matlab/ols/sur.m
+++ b/matlab/ols/sur.m
@@ -46,7 +46,7 @@ jsonmodel = jsonmodel.model;
 %% Find parameters and variable names in equations and setup estimation matrices
 M_exo_names_trim = cellstr(M_.exo_names);
 M_param_names_trim = cellstr(M_.param_names);
-[X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, pidxs] = ...
+[X, Y, startdates, enddates, startidxs, residnames, pbeta, vars, pidxs, surconstrainedparams] = ...
     pooled_sur_common(ds, lhs, rhs, lineno, M_exo_names_trim, M_param_names_trim);
 
 if size(X, 2) ~= M_.param_nbr
@@ -82,7 +82,7 @@ X = newX;
 % Estimated Parameters
 oo_.sur.dof = length(maxfp:minlp);
 [q, r] = qr(X, 0);
-xpxi = (r'*r)\eye(M_.param_nbr);
+xpxi = (r'*r)\eye(size(X, 2));
 resid = Y - X * (r\(q'*Y));
 resid = reshape(resid, oo_.sur.dof, length(lhs));
 
@@ -145,6 +145,12 @@ if ~options_.noprint
         sprintf('s^2: %f', oo_.sur.s2), ...
         sprintf('Durbin-Watson: %f', oo_.sur.dw)};
 
+    if ~isempty(surconstrainedparams)
+        afterward = [afterward, ...
+            sprintf('Constrained parameters: %s', ...
+            strjoin(pbeta(surconstrainedparams), ', '))];
+    end
+
     dyn_table('SUR Estimation', preamble, afterward, [vars{:}], ...
         {'Coefficients','t-statistic','Std. Error'}, 4, ...
         [oo_.sur.beta oo_.sur.tstat oo_.sur.stderr]);