
#Chatterjee, Morley, and Singh (CMS) "Estimating Household Consumption Insurance",
#forthcoming, Journal of Applied Econometrics
#Working paper version available at https://ssrn.com/abstract=2933226
#
#If you use this code, please cite the paper
#
#This version of the qmle_CMS_FD.R code: 11 November 2020
#Output file is qmle_CMS_FD_R.txt with no bias correction and qmle_CMS_bbc_FD_R.txt with bias correction
#
# For a new dataset, it is a good idea to check for different starting
# values on line 70 to get global max. Modify assumptions for different shock
# variances over time for different datasets. Modify the state - space form
# for different models than BPP. Set missing observations in dataset to 999.

rm(list=ls())

library(numDeriv)
library(ucminf)
source("trans.R")
source("lik_fcn_FD.R")
source("scr_fcn_FD.R")
source("shocks_FD.R")
source("dgp_FD.R")
options(warn=-1)

skip_estimation =  0 #set to 1 if estimation already done and just conducting estimation of biases given estimates and 0 otherwise
bootstrap =  0 #set to 1 if using semi-parametric bootstrap to bias-adjust estimates and 0 otherwise
bsim =  50 #set maximum number of bootstrap simulations to estimate biases
number_sim_start =  1 #update this after any crash with estimation for bootstrap simulations to the number_sim at which it crashed

# Select data code here. See below for descriptions. 
data_code <- 1

if(data_code == 1){ # Select for whole sample
    N <- 1765
    Ts <- 15
    data_im <- read.table("cohA_CMS_added1978.txt",sep=",")
    
} else if(data_code == 2){ # Select for low education
    N <- 883
    Ts <- 15 
    data_im <- read.table("cohA_low_CMS_added1978.txt",sep=",")
    
} else if(data_code == 3){ # Select for high education
    N <- 882
    Ts <- 15 
    data_im <- read.table("cohA_high_CMS_added1978.txt",sep=",")
    
} else if(data_code == 4){ # Select for age 1 (30-47 years)
    N <- 1413
    Ts <- 15 
    data_im <- read.table("cohA_age1_CMS_added1978.txt",sep=",")
    
} else if(data_code == 5){ # Select for age 2 (48-65 years)
    N <- 708
    Ts <- 15 
    data_im <- read.table("cohA_age2_CMS_added1978.txt",sep=",")
}

TN = N*Ts
data_sample = cbind(data_im[1:TN, 1], data_im[1:TN, 2]) # Residual income is in column 1 and residual consumption is in column 2. 

start <- 2 # Evaluate likelihood in second quarter of levels data, as in Morley, Nelson, and Zivot (2003).

y = data_sample

## MLE ##############################################################

prmtr_in = c(matrix(-2,36,1)) # Initial parameter values

if (skip_estimation==1){

estimates <- read.table("qmle_CMS_FD_R_store.txt")	
prm_fnl = estimates[1:36,1]
sd_fnl = estimates[1:36,2]
xout = estimates[1:36,3]   
fout = estimates[1,4]	
	
} else{

# Optimization routine
model = ucminf(prmtr_in,lik_fcn_FD, gr = NULL, control=list(trace=1), hessian = 1)
# Returns paramter estimates, -LL value, code, gradient, hessian                                

xout = model$par
fout = model$value

# Final parameter values
prm_fnl = trans(xout)

# Use Hessian to find parameter standard errors
hessn0 = -model$hessian
cov0_h = solve(-hessn0,tol = 1e-35)

grad0 = jacobian(scr_fcn_FD, model$par)
ggg = t(grad0)%*%grad0
cov0_op = solve(ggg,tol = 1e-35)

ihess = solve(-hessn0,tol = 1e-35)
cov0_sw = ihess%*%ggg%*%ihess

grdn_fnl = jacobian(trans, model$par)
cov = grdn_fnl%*%cov0_sw%*%t(grdn_fnl)

sd_fnl = sqrt(diag(cov)) # Standard errors of estimated parameters

fout_store = matrix(-fout,length(prm_fnl),1)

results_store = matrix(c(prm_fnl, sd_fnl, xout, fout_store),length(prm_fnl),4)

sink("qmle_CMS_FD_R_store.txt")
print(results_store)
sink()

}

results = matrix(c(prm_fnl, sd_fnl),length(prm_fnl),2)

# Print table of results into 'qmle_CMS_FD_R.txt'
sink("qmle_CMS_FD_R.txt")
cat("Quasi log likelihood value is ", -fout);
print("Parameter, QMLE parameter estimates, HW standard errors")
print(results)
sink()

if (bootstrap==1) {

smoother = 1 #set to 1 if using Kalman smoother to estimate shocks for bootstrap and 0 for Kalman filter

#Estimate empirical shocks
shocks_v = shocks_FD(xout,1)
shocks_u = shocks_FD(xout,2)
shocks_eps = shocks_FD(xout,3)
shocks_eta = shocks_FD(xout,4)

bias_store = matrix(NA,length(prm_fnl),0)
prmtr_index_store = matrix(NA,length(prm_fnl),0)

if (number_sim_start==1){
store = matrix(NA,2*length(prm_fnl)+1,0)
} else if (number_sim_start>1){
store <- read.csv("qmle_bias_bs_FD.csv", header=F,sep=',')
row.names(store)<-NULL
colnames(store)<-NULL
store=t(t(store))	
bias_store <- store[2:37,]
row.names(bias_store)<-NULL
colnames(bias_store)<-NULL
bias_store=t(t(bias_store)) 	
prmtr_index_store <- store[38:73,]
row.names(prmtr_index_store)<-NULL
colnames(prmtr_index_store)<-NULL
prmtr_index_store=t(t(prmtr_index_store))	
}

number_sim=1
while (number_sim <=  bsim){

   if (number_sim >= number_sim_start){  	
   	# simplist way to print within loop
  print(number_sim)
  Sys.sleep(0.01)
  flush.console()
   }

## Draw sample from bootstrap DGP
y <- dgp_FD(prm_fnl)

col_iter <- 1
while (col_iter <=  ncol(y)){

    row_iter <- 1
    while (row_iter <=  nrow(y)){
        if (data_sample[row_iter,col_iter]==999){
            y[row_iter,col_iter] <- 999
        }
    row_iter <- row_iter + 1
    }

col_iter <- col_iter + 1
}

if (number_sim<number_sim_start){
    number_sim <- number_sim + 1
} else {

#Estimation for bootstrap sample
 # Optimization routine
model = ucminf(prmtr_in,lik_fcn_FD, gr = NULL, hessian = 0)
# Returns paramter estimates, -LL value, code, gradient 

# Final parameter values
prm_fnl_bs = trans(model$par)

bias_store = cbind(bias_store, prm_fnl_bs-prm_fnl)

prmtr_index <- matrix(1,36,1)
dbb0 <- matrix(1,36,1)
dbb1 <- matrix(1,36,1)

ncolbias_store=ncol(bias_store)
ncolbias_store1=ncol(bias_store)-1
ncolbias_store2=ncol(bias_store)-2
if (ncolbias_store>2){
bb0 = rowMeans(bias_store)
bb1 = rowMeans(t(t(bias_store[,1:ncolbias_store1])),1)
bb2 = rowMeans(t(t(bias_store[,1:ncolbias_store2])),1)
dbb0 = abs(bb1-bb0)
dbb1 = abs(bb2-bb1)
}

row_no=1:nrow(prmtr_index)

if (sum(prmtr_index[dbb0 < 0.001 & dbb1 < 0.001]) != nrow(prmtr_index)){
prmtr_index[row_no[dbb0 >=  0.001 | dbb1 >=0.001]] <- matrix(0, length(row_no[dbb0 >=  0.001 | dbb1 >=0.001]), 1)
}
if (sum(prmtr_index)>0){
prmtr_index[row_no[prmtr_index == 1 | prmtr_index_store[,ncol(prmtr_index_store)] == 1]] <- matrix(1,length(prmtr_index[row_no[prmtr_index == 1 | prmtr_index_store[,ncol(prmtr_index_store)] == 1]]),1)
}
prmtr_index_store <- cbind(prmtr_index_store,prmtr_index)
row.names(prmtr_index_store)<-NULL
colnames(prmtr_index_store)<-NULL
prmtr_index_store=t(t(prmtr_index_store)) 

if (number_sim>1){
store <- read.csv("qmle_bias_bs_FD.csv", header=F,sep=',')
row.names(store)<-NULL
colnames(store)<-NULL
store=t(t(store)) 
}

bsim_vec=1:bsim

bs_results=t(cbind(bsim_vec[number_sim],t(prm_fnl_bs-prm_fnl)))
bs_results=rbind(bs_results,prmtr_index)
row.names(bs_results)<-NULL
colnames(bs_results)<-NULL
bs_results=t(t(bs_results)) 

store <- cbind(store,bs_results)
row.names(store)<-NULL
colnames(store)<-NULL
write.table(store, file = "qmle_bias_bs_FD.csv",row.names=F, col.names=FALSE, sep=",")

if (sum(prmtr_index)==length(prm_fnl)){
    number_sim <- bsim
}

number_sim <- number_sim + 1
}

}

prm_fnl_bs = prm_fnl-rowMeans(bias_store) #adjusts the estimates for the bootstrap estimate of bias
sd_fnl_bs = apply(bias_store,1,sd)

results = matrix(c(prm_fnl, prm_fnl_bs, sd_fnl, sd_fnl_bs, rowMeans(bias_store)), length(prm_fnl),5)

sink("qmle_CMS_bbc_FD_R.txt")
cat("Quasi log likelihood value is ", -fout);
print("Parameter, QMLE parameter estimates, Bias-corrected QMLE parameter estimates, HW standard errors, bootstrap standard errors, estimated bias")
print(results)
sink()

} # End of bootstrap bias correction


