rm(list=ls())
#set Seed
set.seed(0814)
#set start time
ptm = proc.time()

#installing and loading neccessary packages
install.packages('nloptr')
library(stats)
library(nloptr)

#defining objective function to be minimized via GMM
eval_f = function(x){
  a = x[1]
  c = x[2]
  d = 0.5*(a*c^(3/2)*sqrt(c*(a^2)+4)+(a^2)*(c^2)+2*c) 
  return((T1-(d/(d+1)))^2 + (T2-((c*d + a*c*d*sqrt(d+1))/((d+1)*(c+1+a*c*sqrt(d+1)))))^2)
}

#defining relationship of 'd' parameter to 'a' and 'c'
df = function(x){
  return(0.5*(x[1]*x[2]^(3/2)*sqrt(x[2]*(x[1]^2)+4)+(x[1]^2)*(x[2]^2)+2*x[2]))
}

#setting values of parameters 'a' and 'c' at which simulations would be done
avals = c(rep(0.75, 4), rep(0.25, 4), rep(0.15, 4), rep(0.85,4) , rep(0.35,4), rep(0.65,4))
cvals = c(rep(0.75, 4), rep(0.25, 4), rep(0.85, 4), rep(0.15,4) , rep(0.65,4), rep(0.35,4))
vals1 = c(rep(0.75,6), rep(0.75,6), rep(df(c(0.75,0.75)),6))
vals2 = c(rep(0.25,6), rep(0.25,6), rep(df(c(0.25,0.25)),6))
vals3 = c(rep(0.15,6), rep(0.85,6), rep(df(c(0.15,0.85)),6))
vals4 = c(rep(0.85,6), rep(0.15,6), rep(df(c(0.85,0.15)),6))
vals5 = c(rep(0.35,6), rep(0.65,6), rep(df(c(0.35,0.65)),6))
vals6 = c(rep(0.65,6), rep(0.35,6), rep(df(c(0.65,0.35)),6))

#setting sizes of simulations
len = rep(c(1000,2500,5000,10000),6)

#vectors to store results from simulation
ahatmean = rep(0,length(avals))
ahatmedian = rep(0,length(avals))
ahatmin = rep(0,length(avals))
ahatmax = rep(0,length(avals))
ahatlb = rep(0,length(avals))
ahatub = rep(0,length(avals))
chatmean = rep(0,length(avals))
chatmedian = rep(0,length(avals))
chatmin = rep(0,length(avals))
chatmax = rep(0,length(avals))
chatlb = rep(0,length(avals))
chatub = rep(0,length(avals))
dhatmean = rep(0,length(avals))
dhatmedian = rep(0,length(avals))
dhatmin = rep(0,length(avals))
dhatmax = rep(0,length(avals))
dhatlb = rep(0,length(avals))
dhatub = rep(0,length(avals))


#simulation begins
for(k in 1:length(avals)){
  
  #specify size
  length = 3000+len[k]
  #specify times to repeat
  times = 5000
  
  #current parameter values
  a = avals[k]
  c = cvals[k]
  
  
  ahat = rep(0, times)
  chat = rep(0, times)
  dhat = rep(0, times)
  
  for(i in 1:times){
    
    #Frechet and Levy values
    U = runif(length+1)
    epsilon = -1/log(U)
    V = runif(length)
    S = 0.5/((qnorm(1-(V/2)))^2)
    
    #MAP process
    Z = rep(0,length)
    Z[1] = max(a*c*sqrt(S[1]*epsilon[length+1]), c*epsilon[1])
    
    for(j in 2:length){
      Z[j] = max(a*c*sqrt(S[j]*Z[j-1]), c*epsilon[j])
    }
    
    #deleting non-stationary part
    ZZ = Z[3000:length]
    
    #summation parts of the objective function
    t1 = exp(-1/(ZZ[2:(length(ZZ))]))
    T1 = mean(t1)
    t2 = exp((-1/ZZ[1:(length(ZZ)-1)])+(-1/ZZ[2:(length(ZZ))]))
    T2 = mean(t2)
    
    #optim function to minimize objective function
    optim = optim(par = c(0.5,0.5), fn = eval_f,
                  method = "L-BFGS-B",
                  lower = c(0,0), upper = c(1,1), control = list(maxit = 1000000))
    
    ahat[i] = optim$par[1]
    chat[i] = optim$par[2]
    dhat[i] = df(c(ahat[i],chat[i]))
  }
  
  ahatmean[k] = mean(ahat)
  ahatmedian[k] = median(ahat)
  ahatmin[k] = min(ahat)
  ahatmax[k] = max(ahat)
  ahatlb[k] = sort(ahat)[times*0.025]
  ahatub[k] = sort(ahat)[times*0.975]
  chatmean[k] = mean(chat)
  chatmedian[k] = median(chat)
  chatmin[k] = min(chat)
  chatmax[k] = max(chat)
  chatlb[k] = sort(chat)[times*0.025]
  chatub[k] = sort(chat)[times*0.975]
  dhatmean[k] = mean(dhat)
  dhatmedian[k] = median(dhat)
  dhatmin[k] = min(dhat)
  dhatmax[k] = max(dhat)
  dhatlb[k] = sort(dhat)[times*0.025]
  dhatub[k] = sort(dhat)[times*0.975]
  
  
}

#table of stored values
Table = cbind(ahatmean, ahatmedian, ahatmin, ahatmax, ahatlb, ahatub, chatmean, chatmedian, chatmin, chatmax, chatlb,
              chatub, dhatmean, dhatmedian, dhatmin, dhatmax, dhatlb, dhatub)

#table of stored values with row labels attached
FinTable = rbind(vals1, Table[1:4,],vals2, Table[5:8,], vals3, Table[9:12,], vals4, Table[13:16,], vals5, Table[17:20,], vals6, Table[21:24,])

FinTable

#time spent
proc.time() - ptm