#### LEARNING POINT: VARIANCE OF R-G STILL TOO HIGH. NEED TO ALLOW FOR AVERAGE MATURITY TO CHANGE ####
rm(list=ls())
library(VARext)
library(debtLimits)
hist <- TRUE
start.year <- 1880
if(hist){
load('data/rmg_hist__1880_est.rdata')
}else{
load('data/rmg_est.rdata')
}
# Load the VAR estimates
#### Set up the stochastic process ####
n.pts <- 2 ; n.dirs <- 4 # (2,5) works
cty <- 'USA' ; lags <- 1
disc <- var.disc(l.cty.mle.rest.1$est[[cty]], n.pts, n.dirs, lb=c(-Inf,0) )
## Try converting to multi-period VAR ##
# maty <- 4
# A.multi <- l.cty.mle.rest.1$est[[cty]]$A %^% maty
# a.multi <- ( diag(2) - A.multi ) %*% l.cty.mle.rest.1$est[[cty]]$mu
# Simga.multi <- 1 / maty * Reduce( '+', lapply( (1:maty)-1,
# function(i) (l.cty.mle.rest.1$est[[cty]]$A) %^% i %*%
# l.cty.mle.rest.1$est[[cty]]$Sigma %*% t(l.cty.mle.rest.1$est[[cty]]$A %^%i ) ) )
# disc <- var.disc(list(a=a.multi,A=A.multi,Sigma=Simga.multi,mu=mu.calc(a.multi,A.multi)),
# n.pts, n.dirs, lb=c(-Inf,0) )
disc <- var.disc(l.cty.mle.rest.4$est[[cty]], n.pts, n.dirs, lb=c(-Inf,0) )
# The dscretized VAR
set.seed(42)
m.idx <- markov_sim( 1e6, disc$trans, 0, disc$n.X )[-(1:1e5)]
m.sim <- disc$X[m.idx+1, ]
# Create the simulation
l.var.est <- var.ols( t(m.sim), lags )
# Check thsurp( d, x[1:5], x[-(1:5)][s.idx], TRUE))e VAR discretization
var.compare <- sapply( list( l.var.est, l.cty.mle.rest.1$est[[cty]] ),
function(x) c( a=x$a, mu=x$mu,A=x$A,
Sigma=x$Sigma[lower.tri(x$Sigma,T)]))
var.compare <- rbind( var.compare, cbind( l.var.est$mu, apply(m.sim,2,mean) ) )
var.compare <- cbind( var.compare, apply( var.compare, 1, function(x) abs(diff(x)) ) )
colnames(var.compare) <- c('data', 'disc', 'diff' )
print(var.compare)
# The comparison to thte data
# pdf('~/Dropbox/2017/research/debtLimits/charts/rg_us_disc.pdf')
plot(disc$X, cex=disc$p.lr * 20, pch=16, ylab='Risk-free rate', xlab='Nominal growth rate')
abline( 0, 1, lty=2 )
# dev.off()
# Plot
mar.dft <- par('mar') ; par(mfrow=c(4,4), mar=c(2,2,2,2))
for( i in 1:disc$n.X ) barplot(disc$trans[i,], main=paste0( 'state ', i) )
par(mfrow=c(1,1), mar=mar.dft)
#### Create the surplus process ####
params <- list( R=1+disc$X[,2]/100, G=1+disc$X[,1]/100,
tri=TRUE, trans=disc$trans )
#### Finish the parameter set-up ####
params$surp.sd <- 1 ; params$lambda <- 0 ; # 0 ;
params$phi <- .8 # .8 # .4
params$cont.type <- 'avg' ; params$diff.method <- "ana" ; params$inner.method <- 'all'
params$d.tri <- TRUE # FALSE # TRUE # Triangulare distribution for surplus shocks
An <- 1 / params$R ; Bn <- rep( -1, length(params$R) ) ; Cn <- An ; def <- matrix(0,1,1)
params$it <- 15 ; params$tol <- .1
# d.init <- sol.nonstoch( params ) # rep( 600, disc$n.X )
### Compute the surplus function from the risk-free solution ###
tgt <- target.moms('USA')
### NEED TO UPDATE THE TARGET MOMENTS FOR EACH COUNTRY ###
rg.dta <- if(hist) hist.read(cty,start.year) else rg.read( cty )
interp <- disc.data.interp( disc, filter(rg.dta[,c('gth','rfr')],rep(1/maty,maty), sides = 1)[-(1:(maty-1)),],
lags, v.date = rg.dta$date[-(1:(maty-1))],
mains=c('Nominal growth rate','Risk free rate',
'Interest-growth differential'))
# Create the parameters for the surplus function consistent with no risk of
# default (so solution is guaraneteed).
# params$v.s.coeff <- rf.surp$v.s.coeff
# params$s.shift <- rf.surp$s.shift
# params$surp.sd <- rf.surp$surp.sd
if(hist){
# params$v.s.coeff <- c(309/4, 334/4, 500/4, -3.88, 20)
ave.s <- mean(rg.dta$pb[],na.rm=TRUE )
params$v.s.coeff <- c(60, 150, 350, 0,0)
}else{
params$v.s.coeff <- c(309, 334, 500, -3.88, 22) # c(324.23, 327.4, 750, -2.88, 20.5)
}
# c(324.23, 327.4, 645, -2.88, 20)
params$s.shift <- rep(0,disc$n.X) # c( 0.55, 0.09, -1.8, -1, 1.41, 0.08, -0.34, 0.57, 0.94, -0.63, 0.87 )
params$surp.sd <- .25 # 2.3657
if( hist){
v.surp <- rg.dta$pb
pb.alt <- (rg.dta$revenue - rg.dta$expenditure) / rg.dta$gdp *100
# Ugly join. Should try to improve...
v.surp[is.na(v.surp)] <- pb.alt[is.na(v.surp)]
rf.surp <- rf.fit( params, rg.dta, rf.debt = params$surp( d, x[1:5], x[-(1:5)][s.idx], TRUE))surp( d, x[1:5], x[-(1:5)][s.idx], TRUE))oeff[3],
rf.level=params$v.s.coeff[5], print.level=1,
init = c( params$v.s.coeff[c(1,2,4)], params$s.shift[-1] ),
v.surp=v.surp, d.init=rg.dta$d[1] )
}else{
rf.surp <- rf.fit( params, tgt$dta, rf.debt = params$v.s.coeff[3],
rf.level=params$v.s.coeff[5], print.level=1,
init = c( params$v.s.coeff[c(1,2,4)], params$s.shift[-1] ) )
}
params$v.s.coeff <- rf.surp$v.s.coeff
s.idx.p <- sapply(1:disc$n.X, function(i) sum( interp$s.idx == i ) ) / nrow(interp$s.idx)
params$s.shift <- rf.surp$s.shift
params$s.shift[s.idx.p==0] <- 0
params$surp.sd <- rf.surp$aad # rf.surp$surp.sd
x.max <- 500
params$d.pts <- 25 ; params$x.mult <- 3
# sol <- sol.wrapper(params, cbind(sol.s$p,sol.s$d), plot.on = FALSE )
params$maxit <- 50
# sol.s <- sol.search( params )
sol <- sol.wrapper(params, cbind(rf.surp$p,rf.surp$d), plot.on = FALSE )
plot.z( sol$p, sol$d, params, xlim=c(0,max(sol$p)*1.2), ylim=c(0,max(sol$p)*1.2) )
# for( b3 in params$v.s.coeff[3] + 1:20 * 10 ){
# params$v.s.coeff[3] <- b3
# message('**************************************')
# message('* params$v.s.coeff[3] = ', params$v.s.coeff[3] )
# message('**************************************')
# sol.new <- sol.wrapper(params, cbind(sol$p,sol$d), plot.on = FALSE )
# plot.z( sol$p, sol$d, params, xlim=c(0,max(sol$p)*1.2), ylim=c(0,max(sol$p)*1.2) )
# if(max(abs(sol.new$err)) < 1e-06){
# sol <- sol.new
# }else{
# break()
# }
# }
# params$v.s.coeff[3] <- 840
# sol <- sol.wrapper(params, cbind(rf.surp$p,rf.surp$d), plot.on = FALSE )
# plot.z( sol$p, sol$d, params, xlim=c(0,max(sol$p)*1.2), ylim=c(0,max(sol$p)*1.2) )
# par(mfrow=c(1,1))
par(mfrow=c(1,1))
plot.surp( params, TRUE, c(0,x.max), TRUE ) #, ylim=c(-15,30) )
# points( tgt$dta$cnlb_gdp_lag[-1], rf.surp$v.surp, pch=16, cex=.75, col=interp$s.idx )
points( rg.dta$d, rf.surp$v.surp, pch=16, cex=.75, col=interp$s.idx )
plot( rg.dta$date, v.surp, type='l', lwd=2, col=2)
lines( rg.dta$date, params$s.shift[interp$s.idx], lwd=2 )
abline(h=0)
# points( tgt$dta$cnlb_gdp_lag, sim[,'s'], pch=16, cex=.75, col=interp$s.idx )
# sol <- sol.wrapper(params, cbind(sol$p,sol$d), plot.on = FALSE )
# sol <- sol.wrapper(params, cbind(rf.surp$p,rf.surp$d), plot.on = FALSE )
# sol <- sol.wrapper(params, plot.on = FALSE )
# d.init <- sol.nonstoch(params)
# d.init[d.init==Inf] <- max(d.init[d.init<Inf])
# sol <- sol.wrapper(params, plot.on = FALSE )
plot.z( sol$p, sol$d, params, xlim=c(0,max(sol$p)*1.2), ylim=c(0,max(sol$p)*1.2) )
# sol <- sol.wrapper(params, cbind(sol.s$p,sol.s$d), plot.on = FALSE )
# stop()
# params.hi <- params
# sol.hi <- sol
# for( inc in seq(.005,.5/4,length.out=10)){
# message( "inc = ", inc )
# params.hi$R <- params$R + inc / 100
# sol.hi <- sol.wrapper(params.hi, cbind(sol.hi$p,sol.hi$d), plot.on = FALSE )
# # Equivalent to shift of half a percentage point or so
# }
params.hi <- params
inc <- -.01
params.hi$R <- params$R + inc / 100
sol.hi <- sol.wrapper(params.hi, cbind(sol$p,sol$d), plot.on = FALSE )
stop()
chg <- -diff(l.thresh.1$USA)
elast.hi <- ( sol$d + ( sol$d - sol.hi$d ) * chg / inc )
pdf('~/Dropbox/2017/research/debtLimits/charts/ts_limit.pdf')
plot( tgt$dta$date[-1], sol$d[interp$s.idx] / 4, type='l',
lwd=2, ylim=range( c(sol$d, elast.hi) )/4,
xlab='Date', ylab='Debt limit' )
lines( tgt$dta$date[-1],
elast.hi[interp$s.idx] / 4,
lwd=2, col='blue' )
legend('bottomright', c('LR R-G at 5% critical value', 'LR R-G at 2.5% critical value'),
bty='n', lwd=2, col=c('black','blue') )
dev.off()
pdf('~/Dropbox/2017/research/debtLimits/charts/surp_int.pdf')
plot( tgt$dta$date[-1], sol$params$s.shift[interp$s.idx],
type='l', lwd=2, xlab='Date', ylab='Surplus intercept' )
abline(h=0,lwd=.5)
dev.off()
plot( (params$R - 1)*100, sol$d / 4, xlab='Nominal interest rate', ylab='Debt limit' )
plot( (params$R - params$G)*100, sol$d / 4, xlab='Interest-growth differential',
ylab='Debt limit' )
plot( (params$G - 1)*100, params$s.shift, xlab='Nominal growth rate', ylab='Surplus intercept' )
plot( (params$R - 1)*100, params$s.shift, xlab='Nominal interest rate', ylab='Surplus intercept' )
plot( (params$R - params$G)*100, params$s.shift, xlab='Interest-growth differential', ylab='Surplus intercept' )
plot( params$s.shift, sol$d / 4, xlab='Surplus intercept', ylab='Debt limit' )
save( tgt, sol, interp, params, elast.hi, file='sol_jul2017.rdata' )
### Compute a discretized simulation
n.sim <- 1e6
n.plot <- 150
n.reg <- 1e6
m.sim <- markov_sim( n.sim, disc$trans, which.max(disc$p.lr)-1, length(disc$p.lr))
sim <- disc$X[m.sim+1,]
colnames(sim) <- c('gth', 'rfr')
varnames <- c( 'Growth', 'Int. rate' )
mu.sim <- apply( sim, 2, mean )
l.var.ols.disc <- var.ols( t(sim), 1 )
l.cty.mle.rest.1$est[[cty]]
var.table( list(data=l.cty.mle.rest.1$est[[cty]], disc=l.var.ols.disc),
file=paste0('~/Dropbox/2017/research/debtLimits/tables/', cty, '_var_tab_disc.tex'),
specnames = c( 'Data', 'Simulation' ), varnames = varnames,
caption = 'Estimated restricted VAR based on data, and VAR estimated from discretized approximation' ,
label=paste0('tab:', cty, '_var_tab_disc'), footer=TRUE )
### Create the IRFs ###
n.irf <- 20
wt.fn <- function( v, a, b, c ){
# Solves for constants s,t such that v = t*a + s*b + (1-s-t)*c
# Where a, b, c, and v are vectors
wts <- solve( cbind( a-c, b-c), v-c )
return( c( wts[1], wts[2], 1-sum(wts) ) )
}
lrv <- var_lr_variance( l.cty.mle.rest.1$est[[cty]]$A, l.cty.mle.rest.1$est[[cty]]$Sigma )
lr.cor <- lrv[1,2] / sqrt( prod( diag( lrv ) ) )
# The long-run correlation
v.targ <- list( R=disc$X[1,]-c(0,.1), G=disc$X[1,]-c(.1,0), both=disc$X[1,] - c(lr.cor,.1) )
# v.targ <- list( R=disc$X[9,]-c(0,.1), G=disc$X[9,]-c(.1,0), both=disc$X[9,] - .1 * c(lr.cor,1) )
# Shocks to R, G and both
l.cand.vert <- list( R=c(1,7,16), G=c(1,10,16), both=c(1,9,2) )
# l.cand.vert <- list( R=c(9,2,3), G=c(9,4,5), both=c(9,4,5) )
# The candidate vertices to put weight on
# shk <- list( R=c(0,1), G=c(1,0), both=1*c(lr.cor,1) )
m.Sigma <- l.cty.mle.rest.1$est[[cty]]$Sigma
cor.innov <- ( m.Sigma[1,2] / sqrt( prod( diag( m.Sigma ) ) ) ) * sqrt( m.Sigma[1,1] )
shk <- list( R=c(0,sqrt(m.Sigma[2,2])), G=c(sqrt(m.Sigma[1,1]),0),
both=c( cor.innov, sqrt(m.Sigma[2,2])) )
shk.dist <- lapply( shk, function(e) e / apply(disc$X,2,sum) )
# Need to figure out how to do a shock to the distribution here
p.init <- lapply( shk, function(x) optim( disc$p.lr[-1],
function(p) sum( ( c( 1-sum(p), p ) %*% disc$X - x - disc$p.lr %*% disc$X ) ^ 2 ) )$par )
p.init <- lapply( p.init, function(p) c(1-sum(p),p) )
l.irf <- list()
for( j in 1:3 ){
irf.p <- matrix( 0, nrow=n.irf+1, ncol=nrow(disc$X) )
# irf.p[1,1] <- 1
# verts <- l.cand.vert[[j]]
# ww <- wt.fn( v.targ[[j]], disc$X[verts[1],], disc$X[verts[2],], disc$X[verts[3],] )
irf.p[ 1, ] <- disc$p.lr
# Initialize the state distribution
irf.p[ 2, ] <- p.init[[j]]
# Initialize the shock
# print( ww )
for( i in 2:n.irf){
irf.p[i+1,] <- irf.p[i,] %*% disc$trans
}
# Markov chain evolution
l.irf[[j]] <- irf.p %*% cbind( disc$X * 4, apply(disc$X,1,diff) * 4, sol$d / 4, params$s.shift ) -
rep( 1, n.irf+1 ) %*% ( irf.p[1,] %*% cbind( disc$X * 4, apply(disc$X,1,diff) * 4, sol$d / 4, params$s.shift ) )
# The IRF is an average across states
}
v.col <- c('blue', 'red', 'darkgreen')
v.lty <- c( 1, 2, 3)
v.nm <- c('G','R','RmG','limit','surp')
v.l <- c( 'topright', 'topright', 'bottomright', 'bottomright', 'bottomright' )
for( j in 1:5 ){
X <- 0:n.irf
Y <- do.call( cbind, lapply(l.irf, function(x) x[,j]) )
pdf(paste0('~/Dropbox/2017/research/debtLimits/charts/irf_', cty, '_', v.nm[j], '.pdf' ) )
plot( range(X), range(Y), type='n', xlab='Quarters', ylab='Response')
# for( i in 1:3 ) lines( X, Y[,i], lwd=2, col=v.col[i] )
for( i in 1:2 ) lines( X, Y[,i], lwd=2, col=v.col[i], lty=v.lty[i] )
abline(h=0,lwd=.5)
# legend(v.l[j], c('R shock', 'G shock', 'Correlated shock'), lwd=2, col=v.col, bty='n')
legend(v.l[j], c('R shock', 'G shock'), lwd=2, col=v.col, bty='n')
dev.off()
}
stop()
### TO DO ###
# - Shut down uncertainty
params.0 <- params
params.0$R <- mean(params$R) # rep( mean(params$R), length(params$R) )
params.0$G <- mean(params$G) # rep( mean(params$G), length(params$G) )
params.0$s.shift <- 0 #
params.0$trans <- matrix(1)
params.0$it <- 150
sol.0 <- sol.wrapper( params.0, cbind( 0, max( sol$d ) ), plot.on = TRUE )
plot.z( sol.0$p, sol.0$d, params.0, xlim=c(0,max(sol.0$p)*1.2), ylim=c(0,max(sol.0$p)*1.2) )
stop()
sol.o <- outer.wrapper( sol, params, print_level = 1 )
#### Create the realized shocks ####
mom.series.plot('USA')
v.surp <- surp.fit(sol.o,interp,params,tgt$dta)
sim <- sim_core( c(1,interp$s.idx), sol.o$d.bar, sol.o$d.grid, sol.o$P, sol.o$Q, params,
tgt$dta$cnlb_gdp_lag[1], TRUE, c(0,0,v.surp) )
surp.ee.sim <- mean( sim[-1,'eps'] )
surp.sd.sim <- sd( sim[-1,'eps'] )
# Mean and sd of surplus shock
# Create the fitted surpluses
sol.err <- outer.err( sol.o, params )
plot.sol( sol.o, xlim=c(min(sol$d)-15,max(sol$d)+15) )
# plot.err( sol.err )
stop()
est <- price.diff.min( params, interp, tgt$dta, sol, sol.o, h.0=c(5,-.5), maxit = 100 )
# est <- price.diff.min( est$sol$params, interp, tgt$dta, est$sol, est$sol.o,
# h.0=c(5,-.5), maxit = 100 )
plot.sol(est$sol.o, xlim=c(min(est$sol$d)-15,max(est$sol$d)+15))
sim.est <- sim_core( c(1,interp$s.idx), est$sol.o$d.bar, est$sol.o$d.grid, est$sol.o$P,
est$sol.o$Q, est$sol$params, tgt$dta$cnlb_gdp_lag[1], TRUE, c(0,0,v.surp) )
# est <- price.diff.min( params, interp, tgt$dta, maxit = 5 )
stop()
plot( tgt$dta$pb_gdp, sim[,'s'], pch=16, col=interp$s.idx,
xlab='Data', ylab='Model', main='Primary balance')
abline(0,1,lty=2)
plot( tgt$dta$date, sim[,'d.prime'], type='l', lwd=2, ylim=c(0,max(sol.o$d.bar)) )
lines( tgt$dta$date, tgt$dta$cnlb_gdp_lag, col='blue', lwd=2 )
lines( tgt$dta$date[-1], sim[-1,'d.bar'], col='red', lwd=2 )
plot( tgt$dta$date, tgt$dta$pb_gdp, type='l', lwd=2, ylim=range(c(tgt$dta$pb_gdp,sim[,'s'])) )
lines( tgt$dta$date, sim[,'s'], col='red', lwd=2 )
abline(h=0)
### Create the long debt prices ###
params.lr.1 <- params ; params.lr.1$lambda <- .75 ; params.lr.1$x.mult <- 5
sol.lr.1 <- outer.wrapper( sol, params.lr.1 )
plot.sol(sol.lr.1, xlim=c(min(sol$d)-40,max(sol$d)+10))
ytm <- sol.o$Q ^ - ( 1 - params$lambda )
Q.1 <- q_hat_mat( sol.o$P, sol.o$d.bar, sol.o$Q, sol.o$Q, sol.o$d.grid, params.lr.1$R, params.lr.1$G,
params.lr.1$s.shift, params.lr.1$lambda, params.lr.1$phi, sol.o$e.grid,
params.lr.1$v.s.coeff, params.lr.1$tri, matrix(0), FALSE, params.lr.1$trans, 2 )
ytm.1 <- Q.1 ^ - ( 1 - params.lr.1$lambda )
# One-year debt
params.lr.5 <- params ; params.lr.5$lambda <- .95
# sol.lr.5 <- outer.wrapper( sol, params.lr.5 )
# plot.sol(sol.lr.5)
Q.5 <- q_hat_mat( sol.o$P, sol.o$d.bar, sol.o$Q, sol.o$Q, sol.o$d.grid, params.lr.5$R, params.lr.5$G,
params.lr.5$s.shift, params.lr.5$lambda, params.lr.1$phi, sol.o$e.grid,
params.lr.5$v.s.coeff, params.lr.5$tri, matrix(0), FALSE, params.lr.5$trans, 2 )
ytm.5 <- Q.5 ^ - ( 1 - params.lr.5$lambda )
# Five-year debt
for( i in 1:disc$n.X){
plot( range(sol.o$d.grid), range(c(0,spd[i,][spd[i,]<Inf],
spd.1[i,][spd.1[i,]<Inf],
spd.5[i,][spd.5[i,]<Inf])), type='n',
lwd=2, xlim=sol$d[i]+c(-8,2), main=paste0("i=",i) )
lines( sol.o$d.grid, spd[i,], type='l', lwd=2, col='black' )
lines( sol.o$d.grid, spd.1[i,], type='l', lwd=2, col='red' )
lines( sol.o$d.grid, spd.5[i,], type='l', lwd=2, col='blue' )
abline(v=sol.o$d.grid, lwd=.5)
}
tp.1 <- ytm.1 - ytm
tp.5 <- ytm.5 - ytm
for( i in 1:disc$n.X){
plot( range(sol.o$d.grid), range(c(0, tp.1[i,][tp.1[i,]<Inf],
tp.5[i,][tp.5[i,]<Inf]), na.rm = TRUE ), type='n',
lwd=2, xlim=sol$d[i]+c(-8,2), main=paste0("i=",i) )
lines( sol.o$d.grid, tp.1[i,], type='l', lwd=2, col='red' )
lines( sol.o$d.grid, tp.5[i,], type='l', lwd=2, col='blue' )
abline(v=sol.o$d.grid, lwd=.5)
}
tp.1.sim <- apply( sim, 1, function(x) approx( sol.o$d.grid, tp.1[x['idx'],],
x['d.prime'] )$y )
tp.5.sim <- apply( sim, 1, function(x) approx( sol.o$d.grid, tp.5[x['idx'],],
x['d.prime'] )$y )
tp <- data.frame( date=tgt$dta$date, yr.1=tgt$dta$Int_1y - tgt$dta$Int_3M,
yr.1.sim=100*tp.1.sim, yr.5=tgt$dta$Int_5y - tgt$dta$Int_3M,
yr.5.sim=100*tp.5.sim, yr.5.1=tgt$dta$Int_5y - tgt$dta$Int_1y,
yr.5.1.sim=100*(tp.5.sim-tp.1.sim))
params.l <- params
params.l$v.s.coeff
### THIS ORDERING LOOKS WRONG ###
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.