#######################################################
#
# DISCUSSION METE KILIC SLOWLY UNFOLDING DISASTERS
#
#######################################################
setwd("C:\\Users\\grouss10\\Dropbox\\Research - Discussions\\discussion Fernie March 2020\\")

# "Wrong" model
mu.c          = .0026
sigma.c       = .0020
mu.eta        = .0200

rho.lambda    = .9933
sigma.lambda  = .0083

lambda.H      = .4167
lambda.L      = .0417

p.00          = 1-0.0001266#0.0017
p.11          = 1-0.0208

risk.aversion = 10
delta         = .999

# "True" formulation
mu.z = sigma.lambda^2/(2*rho.lambda)
Beta = rho.lambda/mu.z

# Mean of the intensity process
library(expm)
transit.mat   = matrix(c(p.00,1-p.00,1-p.11,p.11), 2,2, byrow=T)
statio.probas = (transit.mat %^% 1e3)[1,]
1-statio.probas[1]^12

mean.lambda   = lambda.L * statio.probas[1] + lambda.H * statio.probas[2]


# Simulation of the shocks
sim.L = 100*12
nb.sim = 1e4

set.seed(1234)
consumption.shocks  = matrix(rnorm(sim.L * nb.sim), sim.L, nb.sim)
regime.shocks       = matrix(runif(sim.L * nb.sim), sim.L, nb.sim)
intensity.shocks    = matrix(runif(sim.L * nb.sim), sim.L, nb.sim)
true.inten.shocks   = matrix(runif(sim.L * nb.sim), sim.L, nb.sim)

poisson.jump.shocks = matrix(runif(sim.L * nb.sim), sim.L, nb.sim)
gamma.poisson       = matrix(runif(sim.L * nb.sim), sim.L, nb.sim)

# Initialization
Wrong.Gaussian.model  = array(NA, dim = c(sim.L, 4, nb.sim))
True.Gamma.model      = array(NA, dim = c(sim.L, 4, nb.sim))
Sim.Pois.disasters    = array(0, dim = c(sim.L, 2, nb.sim))

Wrong.Gaussian.model[1,,] = c(
  mu.c - mu.eta * mean.lambda,
  mean.lambda,
  0,
  100
) %*% matrix(1,1,nb.sim)
True.Gamma.model[1,,] = Wrong.Gaussian.model[1,,]

# Actual loop
for(k in 1:nb.sim){
  if(k/100==ceiling(k/100)){print(k)}
  for(i in 2:sim.L){
    
    # Wrong model
    #------------
    # regimes
    Wrong.Gaussian.model[i,3,k] = 
      (Wrong.Gaussian.model[i-1,3,k]==0) * (regime.shocks[i,k] < transit.mat[1,2]) +
      (Wrong.Gaussian.model[i-1,3,k]==1) * (regime.shocks[i,k] < transit.mat[2,2])
    True.Gamma.model[i,3,k]     = Wrong.Gaussian.model[i,3,k]
    
    # intensities
    Wrong.Gaussian.model[i,2,k] = (Wrong.Gaussian.model[i,3,k]*lambda.H +
      (1 - Wrong.Gaussian.model[i,3,k])*lambda.L) * (1 - rho.lambda) +
      rho.lambda * Wrong.Gaussian.model[i-1,2,k] +
      sigma.lambda * sqrt(Wrong.Gaussian.model[i-1,2,k]) * 
      qnorm(intensity.shocks[i,k])
    
    sim.Z.pois = qpois(intensity.shocks[i,k], 
                       lambda = Beta * True.Gamma.model[i-1,2,k])
    True.Gamma.model[i,2,k] = 
      (Wrong.Gaussian.model[i,3,k]*lambda.H +
         (1 - Wrong.Gaussian.model[i,3,k])*lambda.L) * (1 - rho.lambda) +
      qgamma(true.inten.shocks[i,k], shape = sim.Z.pois, scale = mu.z) 
    
    # consumption
    Sim.Pois.disasters[i,1,k] = Sim.Pois.disasters[i-1,1,k]+ 
      qpois(poisson.jump.shocks[i,k],
            lambda = Wrong.Gaussian.model[i-1,2,k])
    Sim.Pois.disasters[i,2,k] = Sim.Pois.disasters[i-1,2,k]+
      qpois(poisson.jump.shocks[i,k],
            lambda = True.Gamma.model[i-1,2,k])
    
    Wrong.Gaussian.model[i,1,k] = mu.c + sigma.c * consumption.shocks[i,k] -
      qgamma(gamma.poisson[i,k], 
             shape = Sim.Pois.disasters[i,1,k] - Sim.Pois.disasters[i-1,1,k], 
             scale = mu.eta)
    
    True.Gamma.model[i,1,k] = mu.c + sigma.c * consumption.shocks[i,k] -
      qgamma(gamma.poisson[i,k], 
             shape = Sim.Pois.disasters[i,2,k]-Sim.Pois.disasters[i-1,2,k], 
             scale = mu.eta)
    
    # consumption level
    Wrong.Gaussian.model[i,4,k] = Wrong.Gaussian.model[i-1,4,k] * 
      exp(Wrong.Gaussian.model[i,1,k])
    True.Gamma.model[i,4,k] = True.Gamma.model[i-1,4,k] * 
      exp(True.Gamma.model[i,1,k])
  }
}


######################################################################
######################################################################
# PLOTS
######################################################################
######################################################################
index = c(28,2)

Wrong.Gaussian.model = Wrong.Gaussian.model[1:1068,,]
True.Gamma.model = True.Gamma.model[1:1068,,]
Sim.Pois.disasters = Sim.Pois.disasters[1:1068,,]
sim.L = 1068

pdf(file = "simulated_paths_unique.pdf", width = 9, height = 9*4/5)
par(mfcol=c(4,2), mai = c(.3, .6,.4,.1))
for(i in index){
  wrong.disasters = NULL
  true.disasters  = NULL
  
  for(time in 1:(sim.L-1)){
    Wrong.trough = time + which.min(log(Wrong.Gaussian.model[(time+1):sim.L,4,i]))
    True.trough  = time + which.min(log(True.Gamma.model[(time+1):sim.L,4,i]))
    
    if(
      (log(Wrong.Gaussian.model[Wrong.trough,4,i])-
          log(Wrong.Gaussian.model[time,4,i]) < -.1) &
      (which.max(log(Wrong.Gaussian.model[time:Wrong.trough,4,i]))==1)
    ){
      wrong.disasters = rbind(wrong.disasters, 
                              c(time, Wrong.trough))
    }
    
    if(
      (log(True.Gamma.model[True.trough,4,i])-
       log(True.Gamma.model[time,4,i]) < -.1) &
      (which.max(log(True.Gamma.model[time:True.trough,4,i]))==1)
    ){
      true.disasters = rbind(true.disasters, 
                              c(time, True.trough))
    }
  }
  ##
  wrong.PTT = NULL
  for(value in unique(wrong.disasters[,2])){
    wrong.PTT = rbind(wrong.PTT,
                      c(min(wrong.disasters[wrong.disasters[,2]==value,1]),
                  value))
  }
  
  ##
  true.PTT = NULL
  for(value in unique(true.disasters[,2])){
    true.PTT = rbind(true.PTT,
                      c(min(true.disasters[true.disasters[,2]==value,1]),
                        value))
  }
  
  plot(Wrong.Gaussian.model[,1,i]*100, type='l', main = "consumption growth ",
       xlab="months", ylab = "monthly % change", 
       ylim = 100*c(
         min(Wrong.Gaussian.model[,1,i]),#True.Gamma.model[,1,i])*100,
         max(Wrong.Gaussian.model[,1,i])#,True.Gamma.model[,1,i])*100
       ), panel.first=grid())
  abline(v=which(Wrong.Gaussian.model[,3,i]==1), col="grey70")
  lines(Wrong.Gaussian.model[,1,i]*100)
  #lines(True.Gamma.model[,1,i]*100, col='red')
  points(c(wrong.PTT), Wrong.Gaussian.model[wrong.PTT,1,i]*100, pch=21, bg="red")
  #points(c(true.PTT), True.Gamma.model[true.PTT,1,i]*100, pch=21, bg="red")
  
  
  plot(log(Wrong.Gaussian.model[,4,i]/100), type='l', main = "consumption level (log)",
       xlab="months", ylab = "", 
       ylim = log(c(
         min(Wrong.Gaussian.model[,4,i]/100),#True.Gamma.model[,4,i])/100,
         max(Wrong.Gaussian.model[,4,i]/100)#,True.Gamma.model[,4,i])/100
       )), panel.first=grid())
  abline(v=which(Wrong.Gaussian.model[,3,i]==1), col="grey70")
  lines(log(Wrong.Gaussian.model[,4,i]/100))
  #lines(log(True.Gamma.model[,4,i]/100), col='red')
  points(c(wrong.PTT), log(Wrong.Gaussian.model[wrong.PTT,4,i]/100), pch=21, bg="red")
  #points(c(true.PTT), log(True.Gamma.model[true.PTT,4,i]/100), pch=21, bg="red")
  
  plot(Sim.Pois.disasters[,1,i], type='l', main = "Consumption jumps",
       xlab="months", ylab = "", 
       ylim = c(
         min(Sim.Pois.disasters[,,i], na.rm = T),
         max(Sim.Pois.disasters[,,i], na.rm = T)
       ), panel.first=grid())
  abline(v=which(Wrong.Gaussian.model[,3,i]==1), col="grey70")
  lines(Sim.Pois.disasters[,1,i])
  #lines(Sim.Pois.disasters[,2,i], col='red')
  #if(i==1){legend("topleft", c("Gaussian", "Gamma"), lty=c(1,1), col=c("black", "red"), bty = "n")}
  
  plot(Wrong.Gaussian.model[,2,i], type='l', main = "Intensities",
       xlab="months", ylab = "", 
       ylim = c(
         min(Wrong.Gaussian.model[,2,i]),#True.Gamma.model[,2,i]),
         max(Wrong.Gaussian.model[,2,i])#,True.Gamma.model[,2,i])
       ), panel.first=grid())
  abline(v=which(Wrong.Gaussian.model[,3,i]==1), col="grey70")
  lines(Wrong.Gaussian.model[,2,i])
  #lines(True.Gamma.model[,2,i], col='red')
}
dev.off()

# Comparing with real data
Quarterly.sim.consgrowth = array(NA, c(sim.L/4, 2, nb.sim))
for(k in 1:nb.sim){
  for(i in 1:(sim.L/4)){
    Quarterly.sim.consgrowth[i,1,k] = sum(Wrong.Gaussian.model[(4*(i-1)+1):(4*i),1,k])
    Quarterly.sim.consgrowth[i,2,k] = sum(True.Gamma.model[(4*(i-1)+1):(4*i),1,k])
  }
}


par(mfrow=c(1,1))
plot(100*Quarterly.sim.consgrowth[,1,1], type='l')
plot(100*Quarterly.sim.consgrowth[,1,6], type='l')
plot(100*log(conso.ratio), type='l')

# Compute dynamics with or without disasters
is.disaster = matrix(0, nb.sim, 2)
for(k in 1:nb.sim){
  
  if(k/100==ceiling(k/100)){print(k)}
  
  for(time in 1:(sim.L-1)){
    wrong.max.PTT = min(log(Wrong.Gaussian.model[time:sim.L,4,k])) - 
      log(Wrong.Gaussian.model[time,4,k])
    true.max.PTT = min(log(True.Gamma.model[time:sim.L,4,k])) - 
      log(True.Gamma.model[time,4,k])
    
    if(wrong.max.PTT < -.1){is.disaster[k,1]=1}
    if(true.max.PTT  < -.1){is.disaster[k,2]=1}
  }

}

# Compute distributions
library(stringr)
dates.real.data = str_c(c(1947:2013)%x%(rep(1,4)),"-", rep(c(1:4),67))
library(zoo)


pdf(file="conso_growth_correc.pdf", height = 3.5, width = 10)
par(mfrow=c(1,3))
plot(as.yearqtr(dates.real.data)[2:length(dates.real.data)],
     100*log(conso.ratio), type='l', main = "Consumption Growth (data)",
     ylab = "quarter % change", xlab = "", 
     ylim = 100*c(
       min(log(conso.ratio), Quarterly.sim.consgrowth[,1,index]),
       max(log(conso.ratio), Quarterly.sim.consgrowth[,1,index])
     ), panel.first= grid())
abline(h=mean(100*log(conso.ratio)), col='red')
abline(h=mean(100*log(conso.ratio))+2*sd(100*log(conso.ratio)), col='red', lty=2)
abline(h=mean(100*log(conso.ratio))-2*sd(100*log(conso.ratio)), col='red', lty=2)

plot(100*Quarterly.sim.consgrowth[1:267,1,index[1]], 
     type='l', main = "simulated path (1)",
     ylab = "quarter % change", xlab = "quarter", 
     ylim = 100*c(
       min(log(conso.ratio), Quarterly.sim.consgrowth[,1,index]),
       max(log(conso.ratio), Quarterly.sim.consgrowth[,1,index])
     ), panel.first= grid())
abline(h=mean(100*Quarterly.sim.consgrowth[1:267,1,index[1]]), col = 'red')
abline(h=mean(100*Quarterly.sim.consgrowth[1:267,1,index[1]])+2*sd(100*Quarterly.sim.consgrowth[1:267,1,index[1]])
       , lty = 2, col = 'red')
abline(h=mean(100*Quarterly.sim.consgrowth[1:267,1,index[1]])-2*sd(100*Quarterly.sim.consgrowth[1:267,1,index[1]])
       , lty = 2, col = 'red')

plot(100*Quarterly.sim.consgrowth[1:267,1,index[2]], 
     type='l', main = "simulated path (2)",
     ylab = "quarter % change", xlab = "quarter", 
     ylim = 100*c(
       min(log(conso.ratio), Quarterly.sim.consgrowth[,1,index]),
       max(log(conso.ratio), Quarterly.sim.consgrowth[,1,index])
     ), panel.first= grid())
abline(h=mean(100*Quarterly.sim.consgrowth[1:267,1,index[2]]), col = 'red')
abline(h=mean(100*Quarterly.sim.consgrowth[1:267,1,index[2]])+2*sd(100*Quarterly.sim.consgrowth[1:267,1,index[2]])
       , lty = 2, col = 'red')
abline(h=mean(100*Quarterly.sim.consgrowth[1:267,1,index[2]])-2*sd(100*Quarterly.sim.consgrowth[1:267,1,index[2]])
       , lty = 2, col = 'red')
dev.off()

pdf(file = "QQplots_correc.pdf", width = 10, height = 3.5)
par(mfrow = c(1,2))
qqplot(100*log(conso.ratio), 100*Quarterly.sim.consgrowth[1:267,1,index[1]],
       main = "QQplot", ylab="simul (% QoQ)", xlab="observed (% QoQ)",
       pch = 21, bg="red", panel.first = grid())
abline(a=0, b=1)
qqplot(100*log(conso.ratio), 100*Quarterly.sim.consgrowth[1:267,1,index[2]],
       main = "QQplot", ylab="simul (% QoQ)", xlab="observed (% QoQ)",
       pch = 21, bg="red", panel.first = grid())
abline(a=0, b=1)
dev.off()

### Table moments
#----------------
all.dynamics.moments = 
  apply(100*Quarterly.sim.consgrowth[,1,], 2, function(x){c(mean(x), sd(x), quantile(x, probs = c(0.05,0.5,0.95)),
                                                    acf(x, plot = F)$acf[c(1,4)+1], acf(x^2, plot = F)$acf[c(1,4)+1])})
no.disasters.moments = 
  apply(100*Quarterly.sim.consgrowth[,1,is.disaster[,1]==0], 2, function(x){c(mean(x), sd(x), quantile(x, probs = c(0.05,0.5,0.95)),
                                                    acf(x, plot = F)$acf[c(1,4)+1], acf(x^2, plot = F)$acf[c(1,4)+1])})

cross.sec.all.dyna.moments = data.frame(apply(all.dynamics.moments, 1, function(x){c(quantile(x,probs = c(.05,.5,.95)))}))
cross.sec.no.disas.moments = data.frame(apply(no.disasters.moments, 1, function(x){c(quantile(x,probs = c(.05,.5,.95)))}))

names(cross.sec.all.dyna.moments) = c("mean", "sd", "5%", "50%", "95%", "AC(1)", "AC(4)", "AC^2(1)", "AC^2(4)")
names(cross.sec.no.disas.moments) = c("mean", "sd", "5%", "50%", "95%", "AC(1)", "AC(4)", "AC^2(1)", "AC^2(4)")

all.dynamics.moments.BP = data.frame(t(all.dynamics.moments))
no.disasters.moments.BP = data.frame(t(no.disasters.moments))
names(all.dynamics.moments.BP) = c("mean", "sd", "5%", "50%", "95%", "AC(1)", "AC(4)", "AC^2(1)", "AC^2(4)")
names(no.disasters.moments.BP) = c("mean", "sd", "5%", "50%", "95%", "AC(1)", "AC(4)", "AC^2(1)", "AC^2(4)")

conso.moments = 
  c(mean(100*log(conso.ratio)), sd(100*log(conso.ratio)), quantile(100*log(conso.ratio), probs = c(0.05,0.5,0.95)),
    acf(100*log(conso.ratio), plot = F)$acf[c(1,4)+1], acf((100*log(conso.ratio))^2, plot = F)$acf[c(1,4)+1])

library(xtable)
xtable(rbind(conso.moments, cross.sec.no.disas.moments,
             cross.sec.all.dyna.moments), digits = c(1,rep(2,5), rep(3,4)))

par(mfrow=c(1,2))
boxplot(no.disasters.moments.BP, main = "No Disasters", ylab = "QoQ % change", pch =4, outline = F, col = "lightblue")
points(c(1:9), conso.moments, col="red", bg = "red", pch = 21)

boxplot(all.dynamics.moments.BP, main = "All samples", ylab = "QoQ % change", pch =4, outline = F, col = "lightblue")
points(c(1:9), conso.moments, col="red", bg = "red", pch = 21)
# PLOT

apply(is.disaster, 2, mean)

pdf(file = "kernel_correc.pdf", width = 10, height = 5)
par(mfrow=c(1,2))

plot(density(100*log(conso.ratio)), main = "kernel densities (no disasters)", xlab = "quarter % change",
     xlim = c(-3,3), panel.first= grid())
lines(density(100*Quarterly.sim.consgrowth[1:267,1,is.disaster[,1]==0]), col='red')
legend("topleft", c("data","simulations"),
       lty =c(1,1), col=c("black", "red"), bty = "n")

plot(density(100*log(conso.ratio)), main = "kernel densities", xlab = "quarter % change",
     xlim = c(-3,3), panel.first= grid())
lines(density(100*Quarterly.sim.consgrowth[1:267,1,]), col='red')
legend("topleft", c("data","simulations"),
       lty =c(1,1), col=c("black", "red"), bty = "n")
dev.off()