# Efficient Implementation of Loop for Bootstrapping Statistics

**URL:** https://discourse.julialang.org/t/efficient-implementation-of-loop-for-bootstrapping-statistics/25675
**Category:** Performance
**Tags:** loops
**Created:** [June 26, 2019, 2:29am UTC](https://discourse.julialang.org/t/efficient-implementation-of-loop-for-bootstrapping-statistics/25675 "2019-06-26T02:29:17Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![slakeon](https://avatars.discourse-cdn.com/v4/letter/s/e274bd/32.png) [@slakeon](https://discourse.julialang.org/u/slakeon)
#### Post date: [June 26, 2019, 2:29am UTC](https://discourse.julialang.org/t/efficient-implementation-of-loop-for-bootstrapping-statistics/25675/1 "2019-06-26T02:29:17Z")

</div>

I am trying to build a loop for a standard recursive-design residual-based nonparametric bootstrap algorithm for a [Vector Autoregression (VAR)](https://en.wikipedia.org/wiki/Vector_autoregression). However, so far, I’m struggling to find an efficient way of writing the algorithm. My current implementation either takes too long to complete a single replication or allocates too much memory.

`npbootstrap` is the main function which generates bootstrap replicates of the slope parameters of the VAR and the lower Cholesky factor of the covariance matrix of the associated residuals. These would be used later on to compute confidence bands for statistics of interest (e.g. impulse response functions).

Below is a toy example with artificial data and parameter estimates.

```julia
using LinearAlgebra, Statistics

pυ = 3 # VAR (lag) order
nυ = 8 # number of variables
Te = 720 # length of vector time series
kζ = nυ * pυ # total number of regressors per equation
Tζ = Te - kζ # degrees of freedom per equation

Y0 = randn(pυ, nυ) # initial conditions
Aυ = randn(nυ * pυ, nυ) # OLS slope parameter estimates
 U = randn(Te, nυ) # residuals

B = 500 # number of bootstrap replications

# main function
function npbootstrap(Y0, Aυ, nυ, pυ, Te, kζ, Tζ, U, B)
    Uc = transform(Te, Tζ, U)
    Ab = Array{Float64,3}(undef, kζ, nυ, B)
    Pb = Array{Float64,3}(undef, nυ, nυ, B)
    for b in 1:B
        Yr = npresample(Y0, Aυ, nυ, pυ, Te, Uc)
        Yr, Zr = regressors(Yr, nυ, pυ, Te)
        Ar = ls(Yr, Zr)
        Ur = resid(Yr, Ar, Zr)
        Σr = covar(Tζ, Ur)
        Pr = cholesky(Symmetric(Σr)).L
        Ab[:, :, b] = Ar
        Pb[:, :, b] = Pr
    end
    Ab, Pb
end

```

```julia
# auxiliary functions
function lag(Yc, nυ, pυ, Te)
    L = Matrix{Float64}(undef, Te, nυ * pυ)
    for i in 1:pυ, j in 1:nυ
        L[:, j + nυ * (i - 1)] = Yc[(pυ + 1 - i):(pυ + Te - i), j]
    end
    L
end

function strim(Yc, pυ)
    Ys = Yc[1:pυ, :]
end

function etrim(Yc, pυ, Te)
    Ye = Yc[(pυ + 1):(pυ + Te), :]
end

function regressors(Yc, nυ, pυ, Te)
    Yb = lag(Yc, nυ, pυ, Te)
     Z = Yb
    Ye = etrim(Yc, pυ, Te)
    Ye, Z
end

function ls(Yc, Z)
    A = Symmetric(Z'Z) \ (Z'Yc)
end

function resid(Yc, A, Z)
    U = Yc - Z * A
end

function covar(Tζ, U)
    Σ = (U'U) / Tζ
end

function transform(Te, Tζ, U)
    Uc = U .- mean(U, dims=1)
    Uc = sqrt(Te / Tζ) * Uc
end

function npresample(Y0, Aυ, nυ, pυ, Te, Uc)
    Yr = Matrix{Float64}(undef, pυ + Te, nυ)
    Yr[1:pυ, :] = Y0
    Ur = Uc[rand(1:Te, Te), :]
    for t in 1:Te
        Yrt = zeros(1, nυ)
        Urt = Ur[t, :]'
        for l in 1:pυ
            Yrtl = Yr[pυ + t - l, :]'
             Aυl = Aυ[(nυ * (l - 1) + 1):(nυ * l), :]
            Yrt += Yrtl * Aυl
        end
        Yr[pυ + t, :] = Yrt + Urt
    end
    Yr
end

```

```julia
using BenchmarkTools

@benchmark npbootstrap(Y0, Aυ, nυ, pυ, Te, kζ, Tζ, U, B)
BenchmarkTools.Trial:
  memory estimate: 1.03 GiB
  allocs estimate: 4777010
  --------------
  minimum time: 658.557 ms (5.69% GC)
  median time: 727.883 ms (5.17% GC)
  mean time: 721.912 ms (6.41% GC)
  maximum time: 767.051 ms (5.34% GC)
  --------------
  samples: 7
  evals/sample: 1

```

```julia
Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = "C:\Users\User1\AppData\Local\atom\app-1.38.2\atom.exe" -a
  JULIA_NUM_THREADS = 4

```

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 26, 2019, 4:40am UTC](https://discourse.julialang.org/t/efficient-implementation-of-loop-for-bootstrapping-statistics/25675/2 "2019-06-26T04:40:42Z")

</div>

Have you looked at [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/) , especially about

1. checking type stability with `@code_warntype`, and

2. profiling your code?

---

<div class="post-metadata">

### Author: ![slakeon](https://avatars.discourse-cdn.com/v4/letter/s/e274bd/32.png) [@slakeon](https://discourse.julialang.org/u/slakeon)
#### Post date: [June 26, 2019, 1:37pm UTC](https://discourse.julialang.org/t/efficient-implementation-of-loop-for-bootstrapping-statistics/25675/3 "2019-06-26T13:37:51Z")

</div>

Sure! I did check some parts of the Performance Tips sections.

By adding the `@views` macro before every loop in the main and the auxiliary functions, memory allocation is halved and computing time is reduced by approx. one third. `@inbounds` improved only marginally on this. Nevertheless, it still stakes too long to compute.

I was wondering if any of the functions can be improved, in particular, regarding either variable preallocation or the order of indexing of the different arrays. Any help would be greatly appreciated.

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 26, 2019, 4:00pm UTC](https://discourse.julialang.org/t/efficient-implementation-of-loop-for-bootstrapping-statistics/25675/4 "2019-06-26T16:00:22Z")

</div>

> [@slakeon](#):
>
> wondering if any of the functions can be improved

This is why I recommended profiling, and checking type stability. You get more useful information if you break up your code into smaller functions.

---

<div class="post-metadata">

### Author: ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)
#### Post date: [June 26, 2019, 4:28pm UTC](https://discourse.julialang.org/t/efficient-implementation-of-loop-for-bootstrapping-statistics/25675/5 "2019-06-26T16:28:00Z")

</div>

Running

```julia
@code_warntype npbootstrap(Y0, Aυ, nυ, pυ, Te, kζ, Tζ, U, B)

```

your type-stability looks OK, so preallocation is the first thing I’d look at. Running 500 bootstrap replicates is allocating a gigabyte of memory, which seems really high. Just skimming your code, I noticed, for instance, the first lines of the `npresample` and `lag` function allocate new arrays and are called at each bootstrap iteration.

As a general suggestion, you might try defining a custom `struct` or two to hold the relevant state variables. I find this can help me reason more easily about what can be preallocated and updated in place…

---

<div class="post-metadata">

### Author: ![slakeon](https://avatars.discourse-cdn.com/v4/letter/s/e274bd/32.png) [@slakeon](https://discourse.julialang.org/u/slakeon)
#### Post date: [June 27, 2019, 11:39pm UTC](https://discourse.julialang.org/t/efficient-implementation-of-loop-for-bootstrapping-statistics/25675/6 "2019-06-27T23:39:19Z")

</div>

I’m going to try to implement that to see if it helps.