---
title: "Abrevaya and Muris - Interval-censoring with fixed effects - Figure 1"
output: html_notebook
---

Setup, including loading the estimator definitions from the R script "01-auxiliary_files.R".

```{r}
library(tidyverse)
library(foreach)
source("01-auxiliary_files.R")
```

Set parameter values.

```{r}
S <- 1000
cp_sim <- cbind(c(2,3),c(60,70))
set.seed('2016')
```

Main loop. 

```{r}
results <- foreach(s = 1:S, .combine = rbind) %:% 
  foreach(n = c(1000), .combine = rbind) %:%
  foreach(sig = c(5),.combine = rbind) %:% 
  foreach(beta  = c(0,1,2,3,4,5,6,7,8,9,10),.combine = rbind) %do% {
    
    if(beta == 0){
      # if(s % 10 == 0){
        print(s)
      # }
    }
    
    # Generate data
    df_s <- icfe_dgp(n=n, beta1 = beta, sigma = sig, cp = cp_sim)
    # Estimate 
    out_log <- icfe(df = df_s, cp = cp_sim)
    
    b_hat <- out_log$beta_hat
    s_hat <- out_log$sigma_hat
    
    b_hat_sum <- 0
    s_hat_sum <- 0
    
    b_OLS <- lm(I(yi2star-yi1star) ~ I(Xi12-Xi11),data=df_s)$coefficient[2]
    
    return(c(s,n,beta,sig,
             b_hat,s_hat,
             b_hat_sum,s_hat_sum,
             b_OLS))
    
  }
```

Transform the results.

```{r}
results <- data.frame(results)
colnames(results) <- c("s","n","beta","sigma",
                       "bhat_logit","sigmahat",
                       "bhat_sum","shat_sum",
                       "b_OLS")
# 
results_summary <- results %>%
  group_by(n,beta,sigma) %>%
  summarise(
    bias_b = mean(bhat_logit) - mean(beta),
    bias_s = mean(sigmahat) - mean(sigma),
    #
    bias_b_sum = mean(bhat_sum) - mean(beta),
    bias_s_sum = mean(shat_sum) - mean(sigma),
    #
    rmse_b = sqrt(mean( (bhat_logit - beta)^2 )),
    rmse_s = sqrt(mean( (sigmahat - sigma)^2 )),
    #
    rmse_b_sum = sqrt(mean( (bhat_sum - beta)^2 )),
    rmse_s_sum = sqrt(mean( (shat_sum - sigma)^2 )),
    #
    rmse_ols = sqrt(mean( (b_OLS - beta)^2 ))
  ) %>%
  mutate(
    eff_rmse = rmse_ols/rmse_b,
    eff_rmse_sum = rmse_ols/rmse_b_sum
  )
```

Show the results.

```{r}
print(results_summary,digits=2)
```

Generate and save the graph.

```{r}
# ggplot(data=results_summary,aes(x=beta,y=rmse_b)) + geom_line() + ylab("RMSE beta-hat")
ggplot(data=results_summary,aes(x=beta,y=eff_rmse)) + geom_line() + 
  ylab("Efficiency") + xlab('beta_0')

#Save it
ggsave("betahat-efficiency.pdf",width = 20, height = 20, units = "cm")
```
