Fit a multivariate regression model to panel
data, assuming different intercepts and common slopes.

Load the sample data.

load('flu')

The dataset array `flu` contains national CDC
flu estimates, and nine separate regional estimates based on Google^{®} query
data.

Extract the response and predictor data.

Y = double(flu(:,2:end-1));
[n,d] = size(Y);
x = flu.WtdILI;

The responses in `Y` are the nine regional
flu estimates. Observations exist for every week over a one-year period,
so *n* = 52. The dimension of the responses corresponds
to the regions, so *d* = 9. The predictors in `x` are
the weekly national flu estimates.

Plot the flu data, grouped by region.

figure;
regions = flu.Properties.VarNames(2:end-1);
plot(x,Y,'x')
legend(regions,'Location','NorthWest')

Fit the multivariate regression model

$${y}_{ij}={\alpha}_{j}+\beta {x}_{ij}+{\epsilon}_{ij},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=1,\dots ,n;\text{\hspace{0.17em}}\text{\hspace{0.17em}}j=1,\dots ,d,$$

with between-region concurrent correlation

$$COV({\epsilon}_{ij},{\epsilon}_{i{j}^{\prime}})={\sigma}_{j{j}^{\prime}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}j=1,\dots ,d.$$

There are *K* = 10 regression coefficients
to estimate: nine intercept terms and a common slope. The input argument `X` should
be an *n*-element cell array of *d*-by-*K* design
matrices.

X = cell(n,1);
for i=1:n
X{i} = [eye(d) repmat(x(i),d,1)];
end
[beta,Sigma] = mvregress(X,Y);

`beta` contains estimates of the *K*-dimensional
coefficient vector

$${\left({\alpha}_{1},{\alpha}_{2},\dots ,{\alpha}_{9},\beta \right)}^{\prime}\text{\hspace{0.17em}}.$$

`Sigma` contains estimates of the *d*-by-*d* variance-covariance
matrix for the between-region concurrent correlations

$$\left(\begin{array}{ccc}{\sigma}_{11}& \dots & {\sigma}_{1,9}\\ \vdots & \ddots & \vdots \\ {\sigma}_{9,1}& \cdots & {\sigma}_{9,9}\end{array}\right)\text{\hspace{0.17em}}.$$

Plot the fitted regression model.

B = [beta(1:d)';repmat(beta(end),1,d)];
xx = linspace(.5,3.5)';
fits = [ones(size(xx)),xx]*B;
figure;
h = plot(x,Y,'x',xx,fits,'-');
for i = 1:d
set(h(d+i),'color',get(h(i),'color'));
end
legend(regions,'Location','NorthWest');

The plot shows that each regression line has a different intercept
but the same slope. Upon visual inspection, some regression lines
appear to fit the data better than others.