$X_t$ is a vector and follows the following Vasicek process. $$ dX_t=(\mu-K\cdot X_t)dt+\Sigma_x\cdot dZ_t \\ $$ what is the variance of $X_t$?
In scalar form the answer is $\frac{\Sigma_x^2}{2\cdot K}\cdot (1-e^{-2Kt})$. But what about in matrix form?
I am reading a Vasicek interest model paper and asked the author for the code to see how he calculated this and below is his matlab code.
[V,D] = eig(K);
sigma_x = chol(Q)';
sigma_y = V \ sigma_x * sigma_x'/(V') ;
sigma = zeros(3,3) ;
for i = 1:3
for j = 1:3
sigma(i,j) = sigma_y(i,j)/(D(i,i) + D(j,j)) *(exp((D(i,i) + D(j,j))/12) -1) ;
end
end
Q = V * sigma * V' ;
I guess I see the general concept but still don't understand why this had to be done element-wise. Also, why isn't there a identity matrix involved, only scalar 1?
Below is my python version of the above code for those who prefer python.
import numpy as np
D, V = np.linalg.eig(K)
D = np.diag(D)
sigma_x = np.linalg.cholesky(Q) # lower triangular matrix
sigma_y = np.linalg.inv(V) @ sigma_x @ sigma_x.T @ np.linalg.inv(V.T)
sigma = np.zeros((3,3))
for i in range(3):
for j in range(3):
sigma[i,j] = sigma_y[i,j] / (D[i,i] + D[j,j]) * (np.exp((D[i,i] + D[j,j])/12) - 1)
Q = V @ sigma @ V.T
12 is there to adjust for time length. Thank you.
Assuming that $K$ and $\Sigma_x$ are square matrices the explicit solution of your vector SDE is $$ X_t=e^{-Kt}\cdot\Big\{X_0+\int_0^t e^{Ks}\cdot\mu\, ds+\Sigma_x\cdot\int_0^t e^{Ks}\cdot dZ_s\Big\}\,. $$ The covariance matrix of the vector $X_t$ is $$ e^{-Kt}\cdot e^{-K^\top t}\cdot\Sigma_x\cdot \Sigma_x^\top\cdot\int_0^te^{Ks}\cdot e^{K^\top s}\,ds\,. $$ When $K$ and $K^\top$ commute and $K+K^\top$ is invertible then this covariance becomes $$ e^{-(K+K^\top)\,t}\cdot\Sigma_x\cdot \Sigma_x^\top\cdot\int_0^te^{(K+K^\top)\,s}\,ds =\Sigma_x\cdot \Sigma_x^\top\cdot(K+K^\top)^{-1}\cdot\Big(1-e^{-(K+K^\top)\,t}\Big)\,. $$