Before proceeding further, we need to address the constraints we placed on the DFA model in Sec 2.2. In particular, we arbitrarily constrained $$\mathbf{Z}$$ in such a way to choose only one of these solutions, but fortunately the different solutions are equivalent, and they can be related to each other by a rotation matrix $$\mathbf{H}$$. Let $$\mathbf{H}$$ be any $$m \times m$$ non-singular matrix. The following are then equivalent DFA models:

$$$\begin{gathered} \mathbf{y}_t = \mathbf{Z}\mathbf{x}_t+\mathbf{a}+\mathbf{v}_t \mathbf{x}_t = \mathbf{x}_{t-1}+\mathbf{w}_t \\ \end{gathered} \tag{10.10}$$$

and

$$$\begin{gathered} \mathbf{y}_t = \mathbf{Z}\mathbf{H}^{-1}\mathbf{x}_t+\mathbf{a}+\mathbf{v}_t \mathbf{H}\mathbf{x}_t = \mathbf{H}\mathbf{x}_{t-1}+\mathbf{H}\mathbf{w}_t \\ \end{gathered}. \tag{10.11}$$$

There are many ways of doing factor rotations, but a common method is the “varimax”" rotation, which seeks a rotation matrix $$\mathbf{H}$$ that creates the largest difference between the loadings in $$\mathbf{Z}$$. For example, imagine that row 3 in our estimated $$\mathbf{Z}$$ matrix was (0.2, 0.2, 0.2). That would mean that green algae were a mixture of equal parts of processes 1, 2, and 3. If instead row 3 was (0.8, 0.1, 0.05), this would make our interpretation of the model fits easier because we could say that green algae followed the first process most closely. The varimax rotation would find the $$\mathbf{H}$$ matrix that makes the rows in $$\mathbf{Z}$$ more like (0.8, 0.1, 0.05) and less like (0.2, 0.2, 0.2).

The varimax rotation is easy to compute because R has a built in function for this: varimax(). Interestingly, the function returns the inverse of $$\mathbf{H}$$, which we need anyway.

## get the estimated ZZ
Z_est <- coef(dfa_1, type = "matrix")$Z ## get the inverse of the rotation matrix H_inv <- varimax(Z_est)$rotmat

We can now rotate both $$\mathbf{Z}$$ and $$\mathbf{x}$$.

## rotate factor loadings
Z_rot = Z_est %*% H_inv
## rotate processes
proc_rot = solve(H_inv) %*% dfa_1\$states