Simple regressions on sliding windows

Simple plots first, to get started.

library(tidyverse)
set.seed(42)
n_sample <- 5000
rho <- 0.95
sd_X <- 0.2
sd_C <- 0.1
sd_r <- 0.3
b <- 1
a <- 0
X <- arima.sim(n = n_sample, list(ar = rho), sd = sd_X) + 1
C <- arima.sim(n = n_sample, list(ar = rho), sd = sd_C)
r <- a + b * C * lead(X) + rnorm(n_sample) * sd_r
data <- data.frame(time = 1:n_sample, X, C, r) 
data[1:100,] %>%
    gather(key = process, value = value, -time) %>%
    ggplot(aes(x = time, y = value, color = process)) + geom_line()

Then, ugly loops. Could be accelerated elegantly via the package slider (see other notebook). The execution is QUITE LONG (~hours).

library(tidyverse)
library(reshape2)

test <- function(t){
    fit <- lm(r ~ C, data = data[(t-j):t,])
    return( (predict(fit, data[t+1,])-data$r[t+1])^2)
}

n_sample <- 10^6
sd_X <- 0.2
sd_C <- 0.1
sd_r <- 0.3
b <- 1
a <- 0
rho <- c(0.8,0.9,0.95,0.99)
wind_max <- 405
wind_min <- 5
wind_by <- 5

err <- array(0, dim= c(length(rho), length(rho), (wind_max-wind_min)/wind_by + 1))
err_inf <- array(0, dim= c(length(rho), length(rho), (wind_max-wind_min)/wind_by + 1))
err_star <- array(0, dim= c(length(rho), length(rho), (wind_max-wind_min)/wind_by + 1))

for(i in 1:length(rho)){
    for(k in 1:length(rho)){
        #if(t%%250==0){print(t)}
        for(j in seq(wind_min, wind_max, by = wind_by)){
            loc <- 1+(j-wind_min)/wind_by
            err_inf[i,k,loc] <- sd_r^2 + sd_X^2 /(1-rho[i]^2) * (a^2 + b^2*sd_C^2/(1-rho[k]^2))
            err_star[i,k,loc] <- sd_r^2 + sd_X^2 * (a^2 + b^2*sd_C^2/(1-rho[k]^2))
        }
    }
}

for(i in 1:length(rho)){
    print(i)
    for(k in 1:length(rho)){
        print(k)
        X <- arima.sim(n = n_sample, list(ar = rho[i]), sd = sd_X) + 1
        C <- arima.sim(n = n_sample, list(ar = rho[k]), sd = sd_C)
        r <- a + b * C * lead(X) + rnorm(n_sample) * sd_r
        data <- data.frame(time = 1:n_sample, X, C, r) 
        for(j in seq(wind_min, wind_max, by = wind_by)){ # T
            loc <- 1+(j-wind_min)/wind_by
            err[i,k,loc] <- map(seq(wind_max+1,n_sample-2, by = 100), test) %>% unlist() %>% mean()
        }
    }
}

err <- melt(err)
err_inf <- melt(err_inf)
err_star <- melt(err_star)
err_full <- cbind(err, err_inf$value, err_star$value)
colnames(err_full) <- c("rho_X", "rho_C", "Window", "Error", "E_inf", "E_star")
for(m in 1:length(rho)){
    err_full$rho_X[err_full$rho_X == m] <- rho[m]
    err_full$rho_C[err_full$rho_C == m] <- rho[m]
}
err_full <- err_full %>%
    mutate(Window = Window * wind_by - 1 + wind_min)
save(err_full, file = "err_full.RData")
# err_full %>% 
#     mutate(rho = paste("rho_X", rho_X, "rho_C", rho_C)) %>%
#     ggplot(aes(x = Window)) + geom_line(aes(y = Error, color = rho)) +
#     xlim(9,100) + ylim(0.08,0.4)

Finally, the plots.

lev <- levels(as.factor(err_full$Window))
lev <- c(lev[1:2], lev[3*(1:(length(lev)/3))])
err_full %>%
    filter(Window %in% lev) %>%
    mutate(rho_X = as.factor(rho_X),
           rho_C = as.factor(rho_C)) %>%
    ggplot(aes(x = Window, y = Error, color = rho_X, shape = rho_C)) +
    geom_point(size = 1.1) + xlim(8,200) + 
    facet_wrap(rho_X ~ rho_C , 
               scales = "free",
               labeller = label_wrap_gen(multi_line=FALSE)) + # + ylim(0.08,0.4)
    geom_line(aes(y = E_inf), color = "black", linetype = 2) + 
    xlab("Window (T)") + ylab("Average Squared Error") +
    scale_color_manual(values = c("#FA7F04", "#2FCD2A", "#0D6DF6", "#F30B58"))
ggsave("err_full.pdf")