%% Implement Seemingly Unrelated Regression
% This example shows how to include exogenous data for several seemingly unrelated regression (SUR) analyses. The response and exogenous series are random paths from a standard Gaussian distribution.
% In seemingly unrelated regression (SUR), each response variable is a function of a subset of the exogenous series, but not of any endogenous variable. That is, for j = 1,..,n and t = 1,...,T, the model for response j at period t is
% y_{jt} = a_{j} + b_{j1}x_{k_1t} + b_{j2}x_{k_2t} + ... + b_{jk_j}x_{k_jt} + \varepsilon_{jt}
% The indices of the regression coefficients and exogenous predictors indicate that:
% You can associate each response with a different subset of exogenous predictors.
% The response series might not share intercepts or regression coefficients.
% SUR accommodates intra-period innovation correlation, but inter-period innovation independence, i.e.,
% E\left( {{\varepsilon _{it}}{\varepsilon _{js}}|X} \right) = \left\{ {\begin{array}{*{20}{c}}
% 0;&t \ne s,\;i\ne j\\
% \sigma _{ij};&i \ne j,\;t = s\\
% {\sigma _i^2 > 0};&i = j,\;t = s
% \end{array}} \right..
%% Simulate Data from True Model
% Suppose that the true model is
% \begin{array}{l}
% {y_{1t}} = 1 + 2{x_{1t}} - 1.5{x_{2t}} + 0.5{x_{3t}} + 0.75{x_{4t}} + \varepsilon_{1t}\\
% {y_{2t}} =  - 1 + 4{x_{1t}} + 2.5{x_{2t}} - 1.75{x_{3t}} - 0.05{x_{4t}} + \varepsilon_{2t}\\
% {y_{3t}} = 0.5 - 2{x_{1t}} + 0.5{x_{2t}} - 1.5{x_{3t}} + 0.7{x_{4t}} + \varepsilon_{3t}
% \end{array},
% where \varepsilon_{jt}, j = 1,...,n are multivariate Gaussian random variables each having mean zero and jointly having covariance matrix
% \Sigma  = \left[ {\begin{array}{*{20}{c}}
% 1&{0.5}&{ - 0.05}\\
% {0.5}&1&{0.25}\\
% { - 0.05}&{0.25}&1
% \end{array}} \right]
% Suppose that the paths represent different econometric measurements, e.g. stock returns.
% Simulate four exogenous predictor paths from the standard Gaussian distribution.

rng(1);   % For reproducibility
n = 3;    % Number of response series 
nExo = 4; % Number of exogenous series
T = 100;
X = randn(100,nExo);

% mvregress, the workhorse of estimate, requires you to input the exogenous data in a T-by-1 cell vector. Cell t of the cell vector is a design matrix indicating the linear relationship of the exogenous variables with each response series at period t. However, estimate associates each predictor to every response. As a result, estimate requires the predictor data in a matrix.
%
% Create a VAR model object that characterizes the true model. Simulate a length 100 path of responses from the model.

aTrue = [1; -1; 0.5];
bTrue = [[2; 4; -2] [-1.5; 2.5; 0.5] [0.5; -1.75; -1.5] [0.75; -0.05; 0.7]];
InnovCov = [1 0.5 -0.05; 0.5 1 0.25; -0.05 0.25 1];
TrueMdl = varm('Beta',bTrue,'Constant',aTrue,'Covariance',InnovCov)

Y = simulate(TrueMdl,T,'X',X);

%
%% SUR Using All Predictors for Each Response Series
% Create a VAR model suitable for SUR using the shorthand syntax of varm.

Mdl1 = varm(n,0);

% Mdl1 is a varm model object template representing a three-dimensional VAR(0) model. Unlike TrueMdl, none of the coefficients, intercepts, and intra-period covariance matrix have values. Therefore, Mdl1 is suitable for estimation.
%
% Estimate the regression coefficients using estimate. Extract the residuals. Display the estimated model using summarize.

[EstMdl1,~,~,E] = estimate(Mdl1,Y,'X',X);
summarize(EstMdl1)

% EstMdl is a varm model object containing the estimated parameters.  E is a T-by- n matrix of residuals.
%
% Alternatively, and in this case, you can use the backslash operator on X and Y. However, you must include a column of ones in X for the intercepts.

coeff = ([ones(T,1) X]\Y)

% coeff is a n-by- nExo + 1 matrix of estimated regression coefficients and intercepts. The estimated intercepts are in the first column, and the rest of the matrix contains the estimated regression coefficients
%
% Compare all estimates to their true values.

InterceptsTbl = table(aTrue,EstMdl1.Constant,coeff(1,:)',...
    'VariableNames',["True" "estimate" "backslash"])
cB = coeff';
cB = cB(:);
CoefficientsTbl = table(bTrue(:),EstMdl1.Beta(:),cB((n + 1):end),...
    'VariableNames',["True" "estimate" "backslash"])

InnovCovTbl = table(InnovCov,EstMdl1.Covariance,...
    'VariableNames',["True" "estimate"])

% The estimates from implementing estimate and the backslash operator are the same, and are fairly close to their corresponding true values.
%
% One way to check the relationship strength between the predictors and responses is to compute the coefficient of determination (i.e., the fraction of variation explained by the predictors), which is
% R^2 = 1 - \frac{\sum_j^n\hat\sigma^2_{\varepsilon j}}{\sum_j^n\hat\sigma^2_{Y j}},
% where \hat\sigma^2_{\varepsilon j} is the estimated variance of residual series j, and \hat\sigma^2_{Y j} is the estimated variance of response series j.

R2 = 1 - sum(diag(cov(E)))/sum(diag(cov(Y)))

% The SUR model and predictor data explain approximately 93% of the variation in the response data.
% Copyright 2012 The MathWorks, Inc.