Packages & parameters

We start by loading the required packages & functions.

rm(list = ls())
library(tidyverse)
library(Rcpp)
library(mvtnorm)
library(reshape2)
library(cowplot)
library(latex2exp)            # For TeX in ggplot
library(viridis)
sourceCpp("c_functions.cpp")
source("R_functions.R")
#folder <- "./sensitivity"

Then, we set some default parameter values. Note: because T is reserved in R, we use TT instead.

Below, we create the grid of parameters that we will test.

a_x <- 0           # Constant in the AR process
a_y <- 0           # Constant in the AR process
n_sim <- 10^6      # Number of simulations

r_x <- c(0.1, 0.5, 0.9, 0.99)
r_y <- c(0.1, 0.5, 0.9)
k   <- c(1, 3, 6, 12, 24)

#r_y <- c(0.1, 0.3, 0.5, 0.7, 0.9, 0.96)
cor <- c(0, 0.5, 0.8)
TT  <- c(6, 12, 24, 36, 60, 84, 120)
pars <- expand.grid(r_x, r_y, cor, k, TT)
r_x <- pars[,1]
r_y <- pars[,2]
cor <- pars[,3]
k <- pars[,4]
TT  <- pars[,5]

Functions & calls

In the tests below, all processes are built so that their variance is equal to one.
The chunk below evaluates the MSE for 8,820 combinations of parameters, it is therefore very time-consuming (hours of CPU), notably because we require one million simulations. Reducing to 100,000 yields more reasonable times. We split the tasks along the k axis.

quad_err <- function(r_x, r_y, cor, TT, k){
    #print(c(r_x, r_y, cor, TT, k))
    quad_err_C_2d2(a_x, a_y, r_x, r_y, 1, 1, cor, TT, k, n_sim)
}

results <- pmap(list(r_x, r_y, cor, TT, k), quad_err) %>% 
    melt() %>% bind_cols(cbind(pars,k))
colnames(results) <- c("MSE", "ID", "r_x", "r_y", "cor", "TT", "k")
save(results, file = "results.RData")

results <- results %>%
    mutate(R2 = 1 - (1-r_y^2) * MSE)

Plot results of sensitivity

Below, we study the sentitivity of the MSE to changes in the parameters set and provide different graphical illustrations.

When k and r_y are independent

First, for k=1.

results %>% 
    filter(k == 1, r_x %in% c(0.1,0.5,0.9,0.99), r_y %in% c(0.1,0.5,0.9), cor %in% c(0,0.5,0.8)) %>%
    ggplot(aes(x = TT, y = R2, color = as.factor(r_x))) + geom_line() +
    facet_grid(cor ~ r_y, scales = "free") +
    # facet_wrap(cor ~ r_y , scales = "free", ncol=3) +
    # theme(strip.text.x = element_blank()) +
    labs(color = TeX("$\\textbf{\\rho_x}$")) +
    xlab(TeX("\\textbf{T}")) + ylab(TeX("\\textbf{R squared}")) +
    theme(legend.position = c(0.2, 0.7)) +
    scale_color_viridis(begin = 0.01, end = 0.91,        # Color management
                        option = "plasma", discrete = T,
                        direction = -1)
ggsave("sensi_k_01.pdf", width = 6, height = 4)    

Then, for k=12.

results %>% 
    filter(k == 12, r_x %in% c(0.1,0.5,0.9,0.99), r_y %in% c(0.1,0.5,0.9), cor %in% c(0,0.5,0.8)) %>%
    ggplot(aes(x = TT, y = R2, color = as.factor(r_x))) + geom_line() +
    facet_grid(cor ~ r_y, scales = "free") +
    # facet_wrap(cor ~ r_y , scales = "free", ncol=3) +
    # theme(strip.text.x = element_blank()) +
    labs(color = TeX("$\\textbf{\\rho_x}$")) +
    xlab(TeX("\\textbf{T}")) + ylab(TeX("\\textbf{R squared}")) +
    theme(legend.position = c(0.92, 0.85)) +
    scale_color_viridis(begin = 0.01, end = 0.91,        # Color management
                        option = "plasma", discrete = T,
                        direction = -1)
ggsave("sensi_k_12.pdf", width = 6, height = 4)    

When k and r_y are linked (memory tradeoff)

We proceeed to a new round of calls.

n_sim <- 10^6
cor <- 0.5 #c(0, 0.5, 0.8)
TT  <- c(1500) # c(6, 12, 24, 36, 60, 84, 120, 150, 200, 250, 300, 500, 1000)
r_x <- 0.9#c(0.1, 0.7, 0.99)
k   <- c(1, 3, 6, 12, 24)#c(1, 2, 3, 4, 5, 6, 12, 24, 36)
pars <- expand.grid(r_x, cor, k, TT)
r_x <- pars[,1]
cor <- pars[,2]
k <- pars[,3]
TT  <- pars[,4]
r_y <- 1-0.9/k
quad_err <- function(r_x, cor, TT, k){
    quad_err_C_2d2(a_x = 0, a_y = 0, 
                   r_x = r_x, r_y = 1 - 0.9/k, 
                   s_x = sqrt(1-r_x^2), s_y = sqrt(1-(1 - 0.9/k)^2), 
                   cor = cor, TT = TT, k = k, n_sim = n_sim)
}

res_sensi <- pmap(list(r_x, cor, TT, k), quad_err) %>% 
    melt() %>% bind_cols(cbind(pars,r_y))
colnames(res_sensi) <- c("value", "ID", "r_x", "cor", "k", "TT", "r_y")

# save(res_sensi, file = "sensi2.RData")


res_sensi <- res_sensi %>%
    mutate(R2 = 1 - (1-r_y^2) * value)

res_sensi %>% 
    filter(TT < 150, cor == 0.5, r_x == 0.9) %>%
    #filter(cor %in% c(0, 0.5, 0.8), r_x %in% c(0.1, 0.7, 0.99), k %in% c(1, 3, 6, 12, 24, 36)) %>%
    ggplot(aes(x = k, y = R2, color = as.factor(TT))) + geom_line() +
    #facet_grid(cor ~ r_x, scales = "free") +
    #facet_wrap(cor ~ r_x, scales = "free", ncol=3) +
    # theme(strip.text.x = element_blank()) +
    labs(color = TeX("$T$")) #+
    # theme(legend.position = c(0.78, 0.84)) 
ggsave("mem.pdf")
Saving 7.29 x 4.51 in image

res_sensi %>% 
    filter(k %in% c(1, 3, 6, 12, 24), TT < 150) %>%
    #filter(cor %in% c(0, 0.5, 0.8), r_x %in% c(0.1, 0.7, 0.99), k %in% c(1, 3, 6, 12, 24, 36)) %>%
    ggplot(aes(x = TT, y = R2, color = as.factor(r_x))) + geom_line() +
    facet_grid(cor ~ r_y, scales = "free") +
    scale_color_viridis(begin = 0.01, end = 0.91,        # Color management
                        discrete = T, option = "plasma", 
                        direction = -1) +
    #facet_wrap(cor ~ r_x, scales = "free", ncol=3) +
    # theme(strip.text.x = element_blank()) +
    labs(color = TeX("$r_x$")) #+
ggsave("sensi.pdf", width = 6, height = 4)

LS0tCnRpdGxlOiAiUHJlZGljdGl2ZSByZWdyZXNzaW9uczogc2Vuc2l0aXZpdHkgYW5hbHlzaXMiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMgUGFja2FnZXMgJiBwYXJhbWV0ZXJzCgpXZSBzdGFydCBieSBsb2FkaW5nIHRoZSByZXF1aXJlZCBwYWNrYWdlcyAmIGZ1bmN0aW9ucy4KCmBgYHtyLCBtZXNzYWdlID0gRiwgd2FybmluZyA9IEZ9CnJtKGxpc3QgPSBscygpKQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShSY3BwKQpsaWJyYXJ5KG12dG5vcm0pCmxpYnJhcnkocmVzaGFwZTIpCmxpYnJhcnkoY293cGxvdCkKbGlicmFyeShsYXRleDJleHApICAgICAgICAgICAgIyBGb3IgVGVYIGluIGdncGxvdApsaWJyYXJ5KHZpcmlkaXMpCnNvdXJjZUNwcCgiY19mdW5jdGlvbnMuY3BwIikKc291cmNlKCJSX2Z1bmN0aW9ucy5SIikKI2ZvbGRlciA8LSAiLi9zZW5zaXRpdml0eSIKYGBgCgpUaGVuLCB3ZSBzZXQgc29tZSAqKmRlZmF1bHQgcGFyYW1ldGVyKiogdmFsdWVzLiAqKk5vdGUqKjogYmVjYXVzZSBUIGlzIHJlc2VydmVkIGluIFIsIHdlIHVzZSBUVCBpbnN0ZWFkLgoKQmVsb3csIHdlIGNyZWF0ZSB0aGUgZ3JpZCBvZiBwYXJhbWV0ZXJzIHRoYXQgd2Ugd2lsbCB0ZXN0LgoKYGBge3J9CmFfeCA8LSAwICAgICAgICAgICAjIENvbnN0YW50IGluIHRoZSBBUiBwcm9jZXNzCmFfeSA8LSAwICAgICAgICAgICAjIENvbnN0YW50IGluIHRoZSBBUiBwcm9jZXNzCm5fc2ltIDwtIDEwXjYgICAgICAjIE51bWJlciBvZiBzaW11bGF0aW9ucwoKcl94IDwtIGMoMC4xLCAwLjUsIDAuOSwgMC45OSkKcl95IDwtIGMoMC4xLCAwLjUsIDAuOSkKayAgIDwtIGMoMSwgMywgNiwgMTIsIDI0KQoKI3JfeSA8LSBjKDAuMSwgMC4zLCAwLjUsIDAuNywgMC45LCAwLjk2KQpjb3IgPC0gYygwLCAwLjUsIDAuOCkKVFQgIDwtIGMoNiwgMTIsIDI0LCAzNiwgNjAsIDg0LCAxMjApCnBhcnMgPC0gZXhwYW5kLmdyaWQocl94LCByX3ksIGNvciwgaywgVFQpCnJfeCA8LSBwYXJzWywxXQpyX3kgPC0gcGFyc1ssMl0KY29yIDwtIHBhcnNbLDNdCmsgPC0gcGFyc1ssNF0KVFQgIDwtIHBhcnNbLDVdCmBgYAoKCiMgRnVuY3Rpb25zICYgY2FsbHMKCkluIHRoZSB0ZXN0cyBiZWxvdywgYWxsIHByb2Nlc3NlcyBhcmUgYnVpbHQgc28gdGhhdCB0aGVpciB2YXJpYW5jZSBpcyBlcXVhbCB0byBvbmUuICAgClRoZSBjaHVuayBiZWxvdyBldmFsdWF0ZXMgdGhlIE1TRSBmb3IgOCw4MjAgY29tYmluYXRpb25zIG9mIHBhcmFtZXRlcnMsIGl0IGlzIHRoZXJlZm9yZSB2ZXJ5IHRpbWUtY29uc3VtaW5nIChob3VycyBvZiBDUFUpLCBub3RhYmx5IGJlY2F1c2Ugd2UgcmVxdWlyZSBvbmUgbWlsbGlvbiBzaW11bGF0aW9ucy4gUmVkdWNpbmcgdG8gMTAwLDAwMCB5aWVsZHMgbW9yZSByZWFzb25hYmxlIHRpbWVzLiAKV2Ugc3BsaXQgdGhlIHRhc2tzIGFsb25nIHRoZSBrIGF4aXMuCgpgYGB7cn0KcXVhZF9lcnIgPC0gZnVuY3Rpb24ocl94LCByX3ksIGNvciwgVFQsIGspewogICAgI3ByaW50KGMocl94LCByX3ksIGNvciwgVFQsIGspKQogICAgcXVhZF9lcnJfQ18yZDIoYV94LCBhX3ksIHJfeCwgcl95LCAxLCAxLCBjb3IsIFRULCBrLCBuX3NpbSkKfQoKcmVzdWx0cyA8LSBwbWFwKGxpc3Qocl94LCByX3ksIGNvciwgVFQsIGspLCBxdWFkX2VycikgJT4lIAogICAgbWVsdCgpICU+JSBiaW5kX2NvbHMoY2JpbmQocGFycyxrKSkKY29sbmFtZXMocmVzdWx0cykgPC0gYygiTVNFIiwgIklEIiwgInJfeCIsICJyX3kiLCAiY29yIiwgIlRUIiwgImsiKQpzYXZlKHJlc3VsdHMsIGZpbGUgPSAicmVzdWx0cy5SRGF0YSIpCgpyZXN1bHRzIDwtIHJlc3VsdHMgJT4lCiAgICBtdXRhdGUoUjIgPSAxIC0gKDEtcl95XjIpICogTVNFKQpgYGAKCiMgUGxvdCByZXN1bHRzIG9mIHNlbnNpdGl2aXR5CgpCZWxvdywgd2Ugc3R1ZHkgdGhlIHNlbnRpdGl2aXR5IG9mIHRoZSBNU0UgdG8gY2hhbmdlcyBpbiB0aGUgcGFyYW1ldGVycyBzZXQgYW5kIHByb3ZpZGUgZGlmZmVyZW50IGdyYXBoaWNhbCBpbGx1c3RyYXRpb25zLgoKCiMjIFdoZW4gayBhbmQgcl95IGFyZSBpbmRlcGVuZGVudAoKRmlyc3QsIGZvciBrPTEuCgpgYGB7cn0KcmVzdWx0cyAlPiUgCiAgICBmaWx0ZXIoayA9PSAxLCByX3ggJWluJSBjKDAuMSwwLjUsMC45LDAuOTkpLCByX3kgJWluJSBjKDAuMSwwLjUsMC45KSwgY29yICVpbiUgYygwLDAuNSwwLjgpKSAlPiUKICAgIGdncGxvdChhZXMoeCA9IFRULCB5ID0gUjIsIGNvbG9yID0gYXMuZmFjdG9yKHJfeCkpKSArIGdlb21fbGluZSgpICsKICAgIGZhY2V0X2dyaWQoY29yIH4gcl95LCBzY2FsZXMgPSAiZnJlZSIpICsKICAgICMgZmFjZXRfd3JhcChjb3IgfiByX3kgLCBzY2FsZXMgPSAiZnJlZSIsIG5jb2w9MykgKwogICAgIyB0aGVtZShzdHJpcC50ZXh0LnggPSBlbGVtZW50X2JsYW5rKCkpICsKICAgIGxhYnMoY29sb3IgPSBUZVgoIiRcXHRleHRiZntcXHJob194fSQiKSkgKwogICAgeGxhYihUZVgoIlxcdGV4dGJme1R9IikpICsgeWxhYihUZVgoIlxcdGV4dGJme1Igc3F1YXJlZH0iKSkgKwogICAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gYygwLjIsIDAuNykpICsKICAgIHNjYWxlX2NvbG9yX3ZpcmlkaXMoYmVnaW4gPSAwLjAxLCBlbmQgPSAwLjkxLCAgICAgICAgIyBDb2xvciBtYW5hZ2VtZW50CiAgICAgICAgICAgICAgICAgICAgICAgIG9wdGlvbiA9ICJwbGFzbWEiLCBkaXNjcmV0ZSA9IFQsCiAgICAgICAgICAgICAgICAgICAgICAgIGRpcmVjdGlvbiA9IC0xKQpnZ3NhdmUoInNlbnNpX2tfMDEucGRmIiwgd2lkdGggPSA2LCBoZWlnaHQgPSA0KSAgICAKYGBgCgpUaGVuLCBmb3Igaz0xMi4KCmBgYHtyfQpyZXN1bHRzICU+JSAKICAgIGZpbHRlcihrID09IDEyLCByX3ggJWluJSBjKDAuMSwwLjUsMC45LDAuOTkpLCByX3kgJWluJSBjKDAuMSwwLjUsMC45KSwgY29yICVpbiUgYygwLDAuNSwwLjgpKSAlPiUKICAgIGdncGxvdChhZXMoeCA9IFRULCB5ID0gUjIsIGNvbG9yID0gYXMuZmFjdG9yKHJfeCkpKSArIGdlb21fbGluZSgpICsKICAgIGZhY2V0X2dyaWQoY29yIH4gcl95LCBzY2FsZXMgPSAiZnJlZSIpICsKICAgICMgZmFjZXRfd3JhcChjb3IgfiByX3kgLCBzY2FsZXMgPSAiZnJlZSIsIG5jb2w9MykgKwogICAgIyB0aGVtZShzdHJpcC50ZXh0LnggPSBlbGVtZW50X2JsYW5rKCkpICsKICAgIGxhYnMoY29sb3IgPSBUZVgoIiRcXHRleHRiZntcXHJob194fSQiKSkgKwogICAgeGxhYihUZVgoIlxcdGV4dGJme1R9IikpICsgeWxhYihUZVgoIlxcdGV4dGJme1Igc3F1YXJlZH0iKSkgKwogICAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gYygwLjkyLCAwLjg1KSkgKwogICAgc2NhbGVfY29sb3JfdmlyaWRpcyhiZWdpbiA9IDAuMDEsIGVuZCA9IDAuOTEsICAgICAgICAjIENvbG9yIG1hbmFnZW1lbnQKICAgICAgICAgICAgICAgICAgICAgICAgb3B0aW9uID0gInBsYXNtYSIsIGRpc2NyZXRlID0gVCwKICAgICAgICAgICAgICAgICAgICAgICAgZGlyZWN0aW9uID0gLTEpCmdnc2F2ZSgic2Vuc2lfa18xMi5wZGYiLCB3aWR0aCA9IDYsIGhlaWdodCA9IDQpICAgIApgYGAKCiMjIFdoZW4gayBhbmQgcl95IGFyZSBsaW5rZWQgKG1lbW9yeSB0cmFkZW9mZikKCldlIHByb2NlZWVkIHRvIGEgbmV3IHJvdW5kIG9mIGNhbGxzLgoKYGBge3IsIGZpZy53aWR0aD00fQpuX3NpbSA8LSAxMF42CmNvciA8LSAwLjUgI2MoMCwgMC41LCAwLjgpClRUICA8LSBjKDE1MDApICMgYyg2LCAxMiwgMjQsIDM2LCA2MCwgODQsIDEyMCwgMTUwLCAyMDAsIDI1MCwgMzAwLCA1MDAsIDEwMDApCnJfeCA8LSAwLjkjYygwLjEsIDAuNywgMC45OSkKayAgIDwtIGMoMSwgMywgNiwgMTIsIDI0KSNjKDEsIDIsIDMsIDQsIDUsIDYsIDEyLCAyNCwgMzYpCnBhcnMgPC0gZXhwYW5kLmdyaWQocl94LCBjb3IsIGssIFRUKQpyX3ggPC0gcGFyc1ssMV0KY29yIDwtIHBhcnNbLDJdCmsgPC0gcGFyc1ssM10KVFQgIDwtIHBhcnNbLDRdCnJfeSA8LSAxLTAuOS9rCnF1YWRfZXJyIDwtIGZ1bmN0aW9uKHJfeCwgY29yLCBUVCwgayl7CiAgICBxdWFkX2Vycl9DXzJkMihhX3ggPSAwLCBhX3kgPSAwLCAKICAgICAgICAgICAgICAgICAgIHJfeCA9IHJfeCwgcl95ID0gMSAtIDAuOS9rLCAKICAgICAgICAgICAgICAgICAgIHNfeCA9IHNxcnQoMS1yX3heMiksIHNfeSA9IHNxcnQoMS0oMSAtIDAuOS9rKV4yKSwgCiAgICAgICAgICAgICAgICAgICBjb3IgPSBjb3IsIFRUID0gVFQsIGsgPSBrLCBuX3NpbSA9IG5fc2ltKQp9CgpyZXNfc2Vuc2kgPC0gcG1hcChsaXN0KHJfeCwgY29yLCBUVCwgayksIHF1YWRfZXJyKSAlPiUgCiAgICBtZWx0KCkgJT4lIGJpbmRfY29scyhjYmluZChwYXJzLHJfeSkpCmNvbG5hbWVzKHJlc19zZW5zaSkgPC0gYygidmFsdWUiLCAiSUQiLCAicl94IiwgImNvciIsICJrIiwgIlRUIiwgInJfeSIpCgojIHNhdmUocmVzX3NlbnNpLCBmaWxlID0gInNlbnNpMi5SRGF0YSIpCgoKcmVzX3NlbnNpIDwtIHJlc19zZW5zaSAlPiUKICAgIG11dGF0ZShSMiA9IDEgLSAoMS1yX3leMikgKiB2YWx1ZSkKYGBgCgoKCmBgYHtyfQoKcmVzX3NlbnNpICU+JSAKICAgIGZpbHRlcihUVCA8IDE1MCwgY29yID09IDAuNSwgcl94ID09IDAuOSkgJT4lCiAgICBnZ3Bsb3QoYWVzKHggPSBrLCB5ID0gUjIsIGNvbG9yID0gYXMuZmFjdG9yKFRUKSkpICsgZ2VvbV9saW5lKCkgKwogICAgbGFicyhjb2xvciA9IFRlWCgiJFQkIikpICMrCmdnc2F2ZSgibWVtLnBkZiIpCgoKcmVzX3NlbnNpICU+JSAKICAgIGZpbHRlcihrICVpbiUgYygxLCAzLCA2LCAxMiwgMjQpLCBUVCA8IDE1MCkgJT4lCiAgICBnZ3Bsb3QoYWVzKHggPSBUVCwgeSA9IFIyLCBjb2xvciA9IGFzLmZhY3RvcihyX3gpKSkgKyBnZW9tX2xpbmUoKSArCiAgICBmYWNldF9ncmlkKGNvciB+IHJfeSwgc2NhbGVzID0gImZyZWUiKSArCiAgICBzY2FsZV9jb2xvcl92aXJpZGlzKGJlZ2luID0gMC4wMSwgZW5kID0gMC45MSwgICAgICAgICMgQ29sb3IgbWFuYWdlbWVudAogICAgICAgICAgICAgICAgICAgICAgICBkaXNjcmV0ZSA9IFQsIG9wdGlvbiA9ICJwbGFzbWEiLCAKICAgICAgICAgICAgICAgICAgICAgICAgZGlyZWN0aW9uID0gLTEpICsKICAgIGxhYnMoY29sb3IgPSBUZVgoIiRyX3gkIikpICMrCmdnc2F2ZSgic2Vuc2kucGRmIiwgd2lkdGggPSA2LCBoZWlnaHQgPSA0KQpgYGAKCg==