---
title: "Abrevaya and Muris - Interval-censoring with fixed effects - Table 2"
output: html_document
---

# Setup

The purpose of this file is to code up a version of the interval-censored fixed effects model with heteroskedasticity, and implement the simulation study that generated the results in Table 2.

```{r}
library(tidyverse)
library(maxLik) # use max lik package
source("01-auxiliary_files.R") # existing estimators

# Some parameters for the design
S <- 1000
cp_sim <- cbind(c(2,3),c(60,70)) # matrix of cut points and correspondence
y_eval <- y_eval2 <- cp_sim[,1]  # cut points
set.seed('2016')                 # RNG seed
```

# Estimator robust to heteroskedasticity

The following modification of the `maxlik` takes into account the heteroskedasticity. 
The starting point for such an analysis is a binary choice model with heteroskedasticity. 
That model is implemented [here](https://www.r-bloggers.com/the-heteroskedastic-probit-model/) and the below simply makes modifications to that code that takes it to the fixed effects case.

```{r}
icfe_CL_het <- function(data,cp_sim,t_start = c(1,1,0)){
  
  # Data is a tibble with variables Xi12, Xi11, yi2, yi1
  # cp_sim is a J-1 by 2 matrix, each row corresponds to
  #    a point in the support of y, except the minimum
  #    of the support.
  # The value of y is in column 1. The corresponding cutoff
  #    is in column 2.
  
  # The difference between this function and the previous one is
  #   that this takes into account the heteroskedasticity in 
  #   the error terms.
  
  # Likelihood function for one cutoff.
  hll_y <- function(cuts,theta){
    
    cut1 <- cuts[1]
    cut2 <- cuts[2]
    
    DX1=(data$Xi12-data$Xi11)
    pX=(data$Xi12+data$Xi11)
    beta <- theta[1]
    gamma0 <- theta[2]
    gamma1 <- theta[3]
    
    d1=ifelse(data$yi1>=y_eval[cut1],1,0)
    d2=ifelse(data$yi2>=y_eval2[cut2],1,0)
    dcut = cp_sim[cut1,2] - cp_sim[cut2,2]
    dbar = d1+d2
    pX = pX[dbar==1]
    DX1 <- DX1[dbar==1]
    d2 <- d2[dbar==1]
    
    q <- 2*d2-1 # sym tranform
    
    het <- exp(gamma0+gamma1*pX)
    
    if(cut1 == cut2){
      xb <- plogis(q*(DX1 * beta)/het)
    } else{
      Z <- cbind(DX1,dcut)  
      xb <- plogis(q*(DX1*beta + dcut)/het)
    }
    
    return(mean(log(xb)))
  }
  
  # Sum acorss likelihoods
  hll <- function(theta){
    
    n_cuts <- nrow(cp_sim)
    X <- expand.grid(1:n_cuts,1:n_cuts)
    return(
      sum(apply(X,1,hll_y,theta = theta))
    )
    
  }
  
  # Maximize composite likelihoods.
  thetahat <- maxLik(hll,start=t_start)$estimate
  
  
  return(
    thetahat
  )
}
```

The following block reproduces one row in Table 2. To reproduce the other rows, change the value of `gamma1`, the parameter governing heteroskedasticity.

```{r, cache=T}
one_sim <- function(s,n = 1000,beta = 1,sig = 2){
  df_s <- icfe_dgp_exp(n=n, beta1 = 1, 
                       gamma0 = log(2), gamma1 = 2,
                       cp = cp_sim)
  
  t1 <- unlist(icfe(df = df_s, cp = cp_sim))[c(2,1)]
  # t2 <- icfe_MD(data = df_s, cp_sim)
  # t3 <- icfe_CL(data = df_s, cp_sim)
  t4 <- icfe_CL_het(data = df_s, cp_sim, t_start = c(t1,0))
  return(
    c(t1,t4)
  ) 
}  
  
out <- as_tibble(t(sapply(1:S,one_sim)))
colnames(out) <- c("CL-b","CL-s","maxlikHet-b","maxlikHet-c0","maxlikHet-c1")

present_out <- out %>% 
  gather(`CL-b`:`maxlikHet-c1`,key="temp",value="estimate") %>% 
  separate(temp, into = c("Estimator","Parameter")) 

# present_out %>% 
  # ggplot() + geom_boxplot(aes(x=Parameter,y=estimate,color=Estimator))

present_out %>% group_by(Parameter,Estimator) %>% 
  summarize(
    bias = mean(estimate),
    rmse = sd(estimate)
  )
```

