library(dlm)
library(doParallel)
cl = makeCluster(detectCores()-1)
registerDoParallel(cl)

t100 <- proc.time()
x <- foreach(j = seq(1,1000),.combine = 'rbind',.packages = c("dlm","doParallel")) %dopar% {

  phi = runif(1,.5,.9)
  sig_epsilon_2 = runif(1,1.1,2)
  print(j)

foreach(i = seq(1,100),.combine='+',.packages = c("dlm","doParallel")) %do% {
        count1_3 = 0
        count2_3 = 0
    print("i")
    print(i)
    print("j")
    print(j)
    t<- proc.time()
    nobs <- 250
    #This is the AR(1) model
    y_t <- arima.sim(n=nobs,list(ar=phi,ma=0),sd=sqrt(sig_epsilon_2))
    #We constrain the 1st parameter to be less than one and
    #the second parameter to be positive.
    parm_rest <- function(parm){
        return( c(exp(parm[1])/(1+exp(parm[1])),exp(parm[2])) )
    }

    original_parameters <- c(phi,sig_epsilon_2)
    ssm_ar1<- function(parm) {
        parm<- parm_rest(parm)
        return(dlm(FF=1,V=0,GG=parm[1],W=parm[2],
        m0=0,C0=solve(1-parm[1]^2)*parm[2]))
    }

    result_3 <- dlmMLE(y_t,parm=c(0,0),build=ssm_ar1,hessian=T)

    coef_3 <- parm_rest(result_3$par)

    dg1_3 <- exp(result_3$par[1])/(1+exp(result_3$par[1]))^2
    dg2_3<- exp(result_3$par[2])
    dg_3<- diag(c(dg1_3,dg2_3))
    myvar_3 <- dg_3%*%solve(result_3$hessian)%*%dg_3
    myvar_3 <- sqrt(diag(myvar_3))

    lower_3 = coef_3 - 1.96 * myvar_3
    upper_3 = coef_3 + 1.96 * myvar_3


############################################################################################################################################
#       count1_3 is the count of the number of times out of 100 that phi falls OUTSIDE it's 95% CI built using dlmMLE
#       count1_3 is ~Binomial(100,.05)  and should have a mean of size * p = 100 * .05 = 5 and
#       a variance of size * p * q = 100 * .05* .95 = 4.75
#       Similarly for count2_3.
#
#       Define Y_n = X_1 + ... + X_n where X_i takes the value 1 in the event that
#       the (generated) original parameter does not lie in it's 95% CI on the ith iteration and 0 otherwise.
#       We compute Y_100 (since i goes from 1 to 100)
#       We can regard this as one realization of Binomial distribution with size=100 and p=.05 (probability of a parameter
#       falling outside it's CI).
#       We repeat the above process j= 1000 times and accumulate the results in count1_1_vector.This will be ~Binomial(n=1000,size=100,p=.05)
#
############################################################################################################################################

        print("Time taken :")
        print(proc.time()-t)

        ifelse(lower_3[1] <= original_parameters[1] && original_parameters[1] <= upper_3[1],print("phi within CI"),
        {print(original_parameters[1]); print("phi outside") ; count1_3= count1_3 +1} )

    ifelse(lower_3[2] <= original_parameters[2] && original_parameters[2] <= upper_3[2],print("sigma_epsilon_2 within CI"),{ print("sigma_epsilon_2 outside") ; count2_3= count2_3 +1} )

    return(c(count1_3,count2_3))


}
}

print(proc.time()-t100)

expected_distribution<- rbinom(1000,100,.05)

count1_3_vector <- x[,1]
count2_3_vector <- x[,2]

pdf("Graph-phi-dlmMLE.pdf",width = 5.6,height = 3.8)
par(mai=c(.8, .8, .3, .2))
qqplot(expected_distribution,count1_3_vector,main=expression(paste("Graph for ",phi)),xlab="Expected distribution",ylab="Observed values")
qqline(count1_3_vector,distribution = function(probs) { qbinom(probs, size=100, prob=0.05) },col = "red",lwd = 2)
dev.off()


pdf("Graph-phi-dlmMLE-jitter.pdf",width = 5.6,height = 3.8)
par(mai=c(.8, .8, .3, .2))
qqplot(jitter(expected_distribution),jitter(count1_3_vector),main=expression(paste("Graph for ",phi," : jittered")),xlab="Expected distribution",ylab="Observed values")
qqline(count1_3_vector,distribution = function(probs) { qbinom(probs, size=100, prob=0.05) },col = "red",lwd = 2)
dev.off()

pdf("Graph-sigma-dlmMLE.pdf",width = 5.6,height = 3.8)
par(mai=c(.8, .8, .3, .2))
qqplot(expected_distribution,count2_3_vector,main=expression(paste("Graph for ",sigma^2)),xlab="Expected distribution",ylab = "Observed values")
qqline(count2_3_vector,distribution = function(probs) { qbinom(probs, size=100, prob=0.05) },col = "red",lwd = 2)
dev.off()


pdf("Graph-sigma-dlmMLE-jitter.pdf",width = 5.6,height = 3.8)
par(mai=c(.8, .8, .3, .2))
qqplot(jitter(expected_distribution),jitter(count2_3_vector),main=expression(paste("Graph for ",sigma^2," : jittered")),xlab="Expected distribution",ylab = "Observed values")
qqline(count2_3_vector,distribution = function(probs) { qbinom(probs, size=100, prob=0.05) },col = "red",lwd = 2)
dev.off()