$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. And 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.
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.
