Generate Covariance Matrix with Pre-Specified Restrictions

Specific Domains Statistics
1

I have a 3 by 3 matrix of random variables and I would like to generate the covariance matrix of these 9 random variables. The restrictions are as follows:

\begin{aligned} \operatorname{Corr}\left(b_{u_{1}, d_{1}}, b_{u_{2}, d_{2}}\right) &=\rho_{1}, \text { if } u_{1} \neq u_{2}, d_{1} \neq d_{2} \\ \operatorname{Corr}\left(b_{u_{1}, d_{1}}, b_{u_{2}, d_{1}}\right) &=\rho_{2}, \text { if } u_{1} \neq u_{2} \\ \operatorname{Corr}\left(b_{u_{1}, d_{1}}, b_{u_{1}, d_{2}}\right) &=\rho_{3}, \text { if } d_{1} \neq d_{2} \\ \operatorname{SD}\left(b_{u_{1}, d_{1}}\right) &=\sigma \end{aligned}

where \rho and \sigma are parameters to be estimated. I write a function to check if two variables are in the same row or column:

# Assign covariances and variances
function fun_cov(p1::Tuple{Int64,Int64}, p2::Tuple{Int64,Int64}, rho1_est::Float64, rho2_est::Float64, rho3_est::Float64, sigma_est::Float64)
    agreements = sum(p1 .== p2)
    ifelse(agreements==2, sigma_est^2, ifelse(agreements==0, rho1_est * sigma_est^2, ifelse(p1[2]==p2[2], rho2_est * sigma_est^2, rho3_est * sigma_est^2)))
end

Then, I generate the covariance matrix:

ind = vec(collect(Iterators.product(1:N, 1:N))) # initialize the index
cov =[fun_cov(p1, p2, rho1_est, rho2_est, rho3_est, sigma_est) for p1 in ind, p2 in ind]

Here is the problem, I find out that some values of \rho and \sigma do not generate a covariance matrix. For example, when

rho1_est=0.6; rho2_est=0.4; rho3_est=0.2; sigma_est=2.0; N=3;

the cov I get has negative eigenvalues. I am not sure what is wrong in my code and I do not see why this setting of correlations is not possible. In my code, I need to guess the value of these parameters to generate the covariance matrix, which is then used to compute the objective function. I wonder how should I make sure cov is indeed a covariance matrix while satisfying the restrictions on correlations. Thanks in advance.

2

Unfortunately it looks like these restrictions don’t specify a valid covariance matrix. You’ll have to change your rho values to specify a positive-definite matrix.

2 Likes
3

Is there any particular reason for this set of restrictions, instead of more or less standard ways of parametrizing covariance matrices (transforming reals into Cholesky factors of correlation/covariance)?

4

thanks, I am thinking about to have a penalty term in my objective function…

5

Those assumptions are implied by the specification of my model.

6

You might want to be careful about your model, then. It’s very easy to write covariance functions that seem reasonable but are not valid. If you happen to have the Brockwell and Davis time series book (Theory and Methods, not the other one), Example 1.5.1 is an example of this that seems close to the issue you’re having: A covariance function for a time series with K(0) = 1, K(\pm 1) = \rho, and everywhere else zero, is only a valid covariance function if |\rho| \leq 1/2. There are some illuminating mathematical arguments about this later in the book, but you could also convince yourself about this by generating the matrix numerically for a sufficiently large data size and using something like the Julia command isposdef.

2 Likes
7

Thanks a lot for the reference. I will look into it. Currently, I am using isposdef to check if it is truly a covariance matrix in my objective function, if not, the objective function returns a large number.