Hi All, thank you so much for all the suggestions, I will try one by one. As many of you suggested, let me post the detail of my code below. So the big picture is an optimization problem called simulated method of moments in which I need to estimate parameters in the multivariate normal distribution. It also involves simulation of this distribution for numerical integration.
I have an N by N matrix of random variables, i.e., there are N^2 random variables in total. Here, I set N=30 for testing the code, but it would be N=500 for the real problem. These N^2 random variables follow a multivariate normal distribution, which is defined by
Corr(\nu_{ij},\nu_{i'j'})=0, \text{ if } i\neq i' \text{ and } j\neq j'
Corr(\nu_{ij},\nu_{i'j'})=\rho_\nu, \text{ if } i=i' \text{ or } j=j' \text{ but not both}
SD(\nu_{ij},\nu_{ij})=\sigma_\nu
where (i,j) is the index of the N by N matrix. So, they have the same standard deviation, and they have the same correlation if the two random variables are in the same column or same row of the matrix. Here are my codes (the issue is in the second function fun_spectral):
using Distributions
using Random
using LinearAlgebra
N = 30 #size of the market
T = 2 #number of markets
num_moments = 2 #number of moments
num_simulation = 2 #number of draws for each numerical integration
rho_nu_est = 0.5 #correlation
sd_nu_est = 1 #standard deviation
# Fix the draws for unobservables.
# I first draw from standard normal distribution, then transform them to multivariate normal.
function fun_fixed_draws(num_moments::Int, num_simulation::Int, T::Int, N::Int)
Random.seed!(1)
#make it an array of arrays, each inner array contain #num_simulation simulations
#the dimension is T * num_moments, so each inner array is for one (T, num_moments) combination
fixed_draws_temp = Array[]
for i = 1:num_moments*T
push!(fixed_draws_temp, randn(N^2, num_simulation))
end
fixed_draws_temp = reshape(fixed_draws_temp, T, num_moments)
return fixed_draws_temp
end
fixed_draws = fun_fixed_draws(num_moments, num_simulation, T, N)
# Transform the draw by spectral method.
# Here, I need to compute the covariance matrix and then factorize it.
function fun_spectral(rho_nu_est, sd_nu_est)
#calculate covariance matrix of $\nu$ given parameters in the distribution
cov_nu_temp = zeros(N^2, N^2)
for i = 1:N^2
cov_nu_temp[i,i] = sd_nu_est^2
end
for i = 1:N
for j = 1:N
for t = 1:N
for k = 1:N
if i == t && j != k
cov_nu_temp[(i-1)*N+j, (t-1)*N+k] = rho_nu_est * sd_nu_est^2
elseif i != t && j == k
cov_nu_temp[(t-1)*N+k, (i-1)*N+j] = rho_nu_est * sd_nu_est^2
end
end
end
end
end
#calculate the spectral decomposition
#I set negative eigenvalues to zero
cov_nu_spectral = eigen(cov_nu_temp)
for i = 1:N^2
if cov_nu_spectral.values[i]<=0
cov_nu_spectral.values[i]=0
end
end
#transform the original draw to draws from multivariate normal distribution
fixed_draws_spec_temp = Array[]
#in the for loop, i first, then t, otherwise the reshape will not put it in correct order
for i = 1:num_moments
for t = 1:T
push!(fixed_draws_spec_temp, cov_nu_spectral.vectors*diagm(sqrt.(cov_nu_spectral.values))*fixed_draws[t,i])
end
end
#reshape draw so that each colum contains the draws for one moments of all markets
fixed_draws_spec_temp = reshape(fixed_draws_spec_temp, T, num_moments)
return fixed_draws_spec_temp
end
fixed_draws_spectral = fun_spectral(rho_nu_est, sd_nu_est)