## 2.9 Solving for the parameters for Form 2*

Solving for the parameters when the model is written in Form 2 is not straight-forward. We could re-write the model in Form 1, or another approach is to use Kronecker products and permutation matrices.

To solve for $$\alpha$$ and $$\beta$$, we need our parameters in a column matrix like so $$\left[ \begin{smallmatrix}\alpha\\\beta\end{smallmatrix} \right]$$. We start by moving the intercept matrix, $$\mathbf{a}$$ into $$\mathbf{Z}$$. $$$\tag{2.24} \begin{bmatrix}stack.loss_1\\stack.loss_2\\stack.loss_3\\stack.loss_4\end{bmatrix} = \begin{bmatrix} \alpha&\beta&0&0&0\\ \alpha&0&\beta&0&0\\ \alpha&0&0&\beta&0\\ \alpha&0&0&0&\beta \end{bmatrix} \begin{bmatrix}1\\air_1\\air_2\\air_3\\air_4\end{bmatrix} + \begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix} = \mathbf{Z}\mathbf{x} + \mathbf{e}.$$$

Then we rewrite $$\mathbf{Z}\mathbf{x}$$ in Equation (2.24) in vec’ form: if $$\mathbf{Z}$$ is a $$n \times m$$ matrix and $$\mathbf{x}$$ is a matrix with 1 column and $$m$$ rows, then $$\mathbf{Z}\mathbf{x} = (\mathbf{x}^\top \otimes \mathbf{I}_n)\,\text{vec}(\mathbf{Z})$$. The symbol $$\otimes$$ means Kronecker product and just ignore it since you’ll never see it again in our course (or google ‘kronecker product’ if you are curious). The “vec” of a matrix is that matrix rearranged as a single column: $\begin{equation*} \,\text{vec} \begin{bmatrix} 1&2\\ 3&4 \end{bmatrix} = \begin{bmatrix} 1\\3\\2\\4 \end{bmatrix} \end{equation*}$ Notice how you just take each column one by one and stack them under each other. In R, the vec is

A = matrix(1:6, nrow = 2, byrow = TRUE)
vecA = matrix(A, ncol = 1)

$$\mathbf{I}_n$$ is a $$n \times n$$ identity matrix, a diagonal matrix with all 0s on the off-diagonals and all 1s on the diagonal. In R, this is simply diag(n).

To show how we solve for $$\alpha$$ and $$\beta$$, let’s use an example with only 3 data points so Equation (2.24) becomes: $$$\tag{2.25} \begin{bmatrix}stack.loss_1\\stack.loss_2\\stack.loss_3\end{bmatrix} = \begin{bmatrix} \alpha&\beta&0&0\\ \alpha&0&\beta&0\\ \alpha&0&0&\beta \end{bmatrix} \begin{bmatrix}1\\air_1\\air_2\\air_3\end{bmatrix} + \begin{bmatrix}e_1\\e_2\\e_3\end{bmatrix}$$$ Using $$\mathbf{Z}\mathbf{x} = (\mathbf{x}^\top \otimes \mathbf{I}_n)\,\text{vec}(\mathbf{Z})$$, this means $$$\begin{bmatrix} \alpha&\beta&0&0\\ \alpha&0&\beta&0\\ \alpha&0&0&\beta \end{bmatrix} \begin{bmatrix}1\\air_1\\air_2\\air_3\end{bmatrix} =\big(\begin{bmatrix}1&air_1&air_2& air_3\end{bmatrix} \otimes \begin{bmatrix}1&0&0\\ 0&1&0 \\ 0&0&1 \end{bmatrix} \bigr) \begin{bmatrix} \alpha\\ \alpha\\ \alpha\\ \beta\\ 0\\ 0\\ 0\\ \beta\\ 0\\ 0\\ 0\\ \beta \end{bmatrix}$$$ We need to rewrite the $$\,\text{vec}(\mathbf{Z})$$ as a permutation’ matrix times $$\left[ \begin{smallmatrix}\alpha\\\beta\end{smallmatrix} \right]$$: $$$\begin{bmatrix} \alpha\\ \alpha\\ \alpha\\ \beta\\ 0\\ 0\\ 0\\ \beta\\ 0\\ 0\\ 0\\ \beta \end{bmatrix} = \begin{bmatrix} 1&0\\ 1&0\\ 1&0\\ 0&1\\ 0&0\\ 0&0\\ 0&0\\ 0&1\\ 0&0\\ 0&0\\ 0&0\\ 0&1\\ \end{bmatrix} \begin{bmatrix} \alpha\\ \beta \end{bmatrix} = \mathbf{P}\mathbf{p}$$$ where $$\mathbf{P}$$ is the permutation matrix and $$\mathbf{p}=\left[ \begin{smallmatrix}\alpha\\\beta\end{smallmatrix} \right]$$. Thus, $$$\tag{2.26} \mathbf{y}=\mathbf{Z}\mathbf{x}+\mathbf{e} = (\mathbf{x}^\top \otimes \mathbf{I}_n)\mathbf{P}\begin{bmatrix}\alpha\\ \beta\end{bmatrix} = \mathbf{M}\mathbf{p} + \mathbf{e}$$$ where $$\mathbf{M}=(\mathbf{x}^\top \otimes \mathbf{I}_n)\mathbf{P}$$. We can solve for $$\mathbf{p}$$, the parameters, using $(\mathbf{M}^\top\mathbf{M})^{-1}\mathbf{M}^\top\mathbf{y}$ as before.

#### 2.9.0.1 Code to solve for parameters in Form 2

In the homework, you will use the R code in this section to solve for the parameters in Form 2.

#make your y and x matrices
y=matrix(dat$stack.loss, ncol=1) x=matrix(c(1,dat$Air.Flow),ncol=1)
#make the Z matrix
n=nrow(dat) #number of rows in our data file
k=1
#Z has n rows and 1 col for intercept, and n cols for the n air data points
#a list matrix allows us to combine "characters" and numbers
Z=matrix(list(0),n,k*n+1)
Z[,1]="alpha"
diag(Z[1:n,1+1:n])="beta"
#this function creates that permutation matrix for you
P=MARSS:::convert.model.mat(Z)$free[,,1] M=kronecker(t(x),diag(n))%*%P solve(t(M)%*%M)%*%t(M)%*%y  [,1] alpha -11.6159170 beta 0.6412918 coef(lm(dat$stack.loss ~ dat$Air.Flow))  (Intercept) dat$Air.Flow
-11.6159170    0.6412918 

Go through this code line by line at the R command line. Look at Z. It is a list matrix that allows you to combine numbers (the 0s) with character string (names of parameters). Notice that class(Z[1,3])="numeric" while class(Z[1,2])="character". This is important. 0 in R is a number while "0" would be a character (the name of a parameter). Look at the permutation matrix P. Try MARSS:::convert.model.mat(Z)$free and see that it returns a 3D matrix, which is why the [,,1] appears (to get us a 2D matrix). To use more data points, you can redefine dat to say dat=stackloss to use all 21 data points. Here’s another example. Rewrite the model with multiple intercepts (Equation (2.14) ) as $$$\tag{2.27} \begin{bmatrix}stack.loss_1\\ stack.loss_2\\ stack.loss_3\\ stack.loss_4\end{bmatrix} = \begin{bmatrix} \alpha_n&\beta&0&0&0\\ \alpha_s&0&\beta&0&0\\ \alpha_n&0&0&\beta&0\\ \alpha_s&0&0&0&\beta \end{bmatrix}\begin{bmatrix}1\\air_1\\air_2\\air_3\\air_4\end{bmatrix} + \begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}=\mathbf{Z}\mathbf{x}+\mathbf{a}+\mathbf{e}$$$ To estimate the parameters, we need to be able to write a list matrix that looks like $$\mathbf{Z}$$ in Equation (2.27). We can use the same code as above with $$\mathbf{Z}$$ changed to look like that in Equation (2.27). y = matrix(dat$stack.loss, ncol = 1)
x = matrix(c(1, dat$Air.Flow), ncol = 1) n = nrow(dat) k = 1 # list matrix allows us to combine numbers and character # strings Z = matrix(list(0), n, k * n + 1) Z[seq(1, n, 2), 1] = "alphanorth" Z[seq(2, n, 2), 1] = "alphasouth" diag(Z[1:n, 1 + 1:n]) = "beta" P = MARSS:::convert.model.mat(Z)$free[, , 1]
M = kronecker(t(x), diag(n)) %*% P
solve(t(M) %*% M) %*% t(M) %*% y
                 [,1]
alphanorth -2.0257880
alphasouth -5.5429799
beta        0.5358166

Similarly to estimate the parameters for Equation (2.18), we change the $$\beta$$’s in our $$\mathbf{Z}$$ list matrix to have owner designations:

Z = matrix(list(0), n, k * n + 1)
Z[seq(1, n, 2), 1] = "alphanorth"
Z[seq(2, n, 2), 1] = "alphasouth"
diag(Z[1:n, 1 + 1:n]) = rep(c("beta.s", "beta.a"), n)[1:n]
solve(t(M) %*% M) %*% t(M) %*% y
            [,1]
beta.a       0.5
The parameters estimates are the same as with the model in Form 1, though $$\beta$$’s are given in reversed order simply due to the way convert.model.mat() is ordering the columns in Form 2’s $$\mathbf{Z}$$.