# Variance-Covariance matrix with missing data

**URL:** <https://discourse.julialang.org/t/variance-covariance-matrix-with-missing-data/83312>\
**Category:** Performance\
**Tags:** statistics, missing-values\
**Created:** [June 24, 2022, 2:10pm UTC](https://discourse.julialang.org/t/variance-covariance-matrix-with-missing-data/83312 "2022-06-24T14:10:41Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Baba\_Yara\_Fahiz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baba_yara_fahiz/32/36478_2.png) [@Baba\_Yara\_Fahiz](https://discourse.julialang.org/u/Baba_Yara_Fahiz)\
**Post date:** [June 24, 2022, 2:10pm UTC](https://discourse.julialang.org/t/variance-covariance-matrix-with-missing-data/83312/1 "2022-06-24T14:10:42Z")

</div>

I am working with financial data where I have 600 matrices that are all 61 by 25000.  
I am computing the variance-covariance matrix but running into performance issues because of missing data.  
The way I am handling the missing data is to compute pairwise covariances across columns.  
Is there any way I can speed this up?  
I have access to a 128 epic system with 220 GB of memory.

```julia
using BenchmarkTools, Random, Missings, StatsBase

M = Matrix{Union{Float64, Missing}}(undef,600,25000)
M .= rand(600,25000) 

ix = rand(CartesianIndices(M), 20_000)  
M[ix] .= missing

function my_cov(x)
    nc = size(x,2);
    t = zeros(nc, nc)
    Threads.@threads for i in 1:nc
        for j in 1:nc
            if i <= j
               sx, sy = skipmissings(x[:, i], x[:, j])
                t[i, j] = cov(collect(sx), collect(sy))
                t[j, i] = t[i, j]
            end 
        end
    end
    return t
end

@benchmark my_cov(M) 

```

Thanks a lot for the suggestions

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [June 24, 2022, 6:01pm UTC](https://discourse.julialang.org/t/variance-covariance-matrix-with-missing-data/83312/2 "2022-06-24T18:01:59Z")

</div>

Just in case: would it be possible to use NaNs instead of missings?  
The function `nancov()` in the [NaNStatistics.jl](https://github.com/brenhinkeller/NaNStatistics.jl) package is very fast.

---

<div class="post-metadata">

**Author:** ![Paul\_Soderlind](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_soderlind/32/1753_2.png) [@Paul\_Soderlind](https://discourse.julialang.org/u/Paul_Soderlind)\
**Post date:** [June 24, 2022, 6:15pm UTC](https://discourse.julialang.org/t/variance-covariance-matrix-with-missing-data/83312/3 "2022-06-24T18:15:11Z")

</div>

Good idea.

In case you want something different, here is an idea:  
(a) for each pair replace `(x,y)` with 0 if either `x` or `y` contains a missing/NaN  
(b) keep track of the number of non-missing/NaN and adjust the covariance estimate accordingly.

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [June 24, 2022, 6:19pm UTC](https://discourse.julialang.org/t/variance-covariance-matrix-with-missing-data/83312/4 "2022-06-24T18:19:40Z")

</div>

I think you want `pairwise` which lives in StatsBase.jl. But I’m not 100% sure how it works. You should read the docs and see if it helps you.

---

<div class="post-metadata">

**Author:** ![Baba\_Yara\_Fahiz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baba_yara_fahiz/32/36478_2.png) [@Baba\_Yara\_Fahiz](https://discourse.julialang.org/u/Baba_Yara_Fahiz)\
**Post date:** [June 24, 2022, 8:20pm UTC](https://discourse.julialang.org/t/variance-covariance-matrix-with-missing-data/83312/5 "2022-06-24T20:20:35Z")

</div>

Thanks a lot for the suggestions.

The pairwise solution is this and runs just as fast as my solution above.

```julia
@benchmark pairwise(cov, eachcol(M), skipmissing =:pairwise, symmetric =true) 

```

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [June 25, 2022, 12:17am UTC](https://discourse.julialang.org/t/variance-covariance-matrix-with-missing-data/83312/6 "2022-06-25T00:17:39Z")

</div>

I know this was not your question, but typically people try not to compute covariance matrices this large as they may not be very trustworthy. You might want to look into shrinkage or other high dimensional estimation techniques. Just an FYI