####################################################### # # 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()