#####################################################################
# plotting.R
#
# Contains the plotting functions
# 22feb2017
# Philip Barrett, Washington DC
#####################################################################
plot.surp <- function(params, non.stoch=TRUE, x.lim=c(0,200), vert=FALSE,
leg=TRUE, ... ){
# Plots the surplus function
d <- seq( from=x.lim[1], to=x.lim[2], by = 1 )
# X values
surps <- sapply( params$s.shift, function(shift)
sapply(d, surp, coeff=params$v.s.coeff, shift=shift, tri=params$tri) )
# Y-values
y.max <- max(surps)
d.max <- d[which.max(if(is.null(dim(surps))) surps else surps[,1])]
if(d.max==0) d.max <- x.lim[2]
y.min <- if( min(surps) < 0 ) ( if( params$tri ) min(surps) else min(surps[d<d.max]) ) else 0
y.lim <- c(y.min,y.max)
# Y range
plot( x.lim, y.lim, type='n', xlab='Debt', ylab='Surplus', ... )
for( i in 1:ncol(surps) ){
lines( d, surps[,i], lwd=2, col=i )
if( non.stoch ){
# sapply( 1:ncol(surps), function(j) abline( 0, params$R[i] - params$G[j], lwd=.5, col=i, lty=2 ) )
abline( 0, params$R[i] - params$G[i], lwd=1, col=i )
}
}
if(leg) legend( 'topright', paste0( 'G = ', round( params$G, 3 ) ), lwd=2, col=1:i, bty='n' )
abline(h=0)
if( vert ){
abline( v=sol.nonstoch(params), lwd=.5, col=1:length(params$R) )
}
}
plot.z <- function( p, d, params, An=c(0), Bn=c(0), Cn=c(0), def=matrix(0), ... ){
# Plots the zed function for all the values of i
n <- length(params$R)
n.x <- ceiling( sqrt(n) )
n.y <- ceiling( n / n.x )
par(mfrow=c(n.y,n.x))
global.apx <- p_init_d( params, p, d, An, Bn, Cn, def )
for( i in 1:n )
plot.z.i( p, d, params, i, An, Bn, Cn, def, global.apx[i], ... )
par(mfrow=c(1,1))
}
plot.z.i <- function( p, d, params, i, An, Bn, Cn, def, global, ylim=NULL,
cex.lab=1, cex.axis=1, cex.main=1, ... ){
# Plots the z function vs. p[i] in state i
p.seq <- c( seq(0,1e-3,by=1e-5), seq(1e-3,1e-2,by=1e-4), seq(1e-2,1e-1,by=1e-3), seq(1e-1,1,by=1e-2) )
# The x-values
y.vals <- t( sapply( p.seq, function(p.i){
this.p <- p
this.p[i] <- p.i
return( zed_2(this.p, d, params, An, Bn, Cn, def )[i,] )
} ) )
# The y values
# if( !exists('ylim')) ylim <- c(0,1)
plot( p.seq, y.vals[,1], lwd=2, xlab='p', ylab='z',
type='l', ylim=ylim, main=paste0( 'i = ', i ),
cex.lab=cex.lab, cex.axis=cex.axis, cex.main=cex.main, ... )
abline( 0, 1, lty=2 )
abline( v=p[i], lty=2 )
abline( h=p[i], lty=2 )
# The level
abline( v=global, lty=2, col='blue' )
# The global minimum
par(new = TRUE)
y.range <- c( max( min(c(0, y.vals[,2]), na.rm = TRUE), -10 ),
min( max(y.vals[,2], na.rm = TRUE ), 10 ) )
# The y limits for the derivative
plot( p.seq, y.vals[,2], type = "l", col='red', axes = FALSE, bty = "n",
xlab = "", ylab = "", ylim=y.range, ... )
# Plot the derivative
axis(side=4, at = pretty(y.range), cex.lab=cex.lab, cex.axis=cex.axis )
# mtext("z.p", side=4, line=3)
abline( h=1, lty=2, col='red' )
# Axis labelling
}
plot.z.d <- function( p, d, params, An=c(0), Bn=c(0), Cn=c(0), def=matrix(0),
d.range=NULL, global=NULL, ... ){
# Plots the zed function for all the values of i
n <- length(params$R)
n.x <- ceiling( sqrt(n) )
n.y <- ceiling( n / n.x )
par(mfrow=c(n.y,n.x))
if(is.null(d.range)) d.range <- c( 0, max(400, max(d)))
# global.apx <- p_init_d( params, p, d, An, Bn, Cn, def )
for( i in 1:n )
plot.z.d.i( p, d, params, i, An, Bn, Cn, def, d.range, global, ... )
par(mfrow=c(1,1))
}
plot.z.d.i <- function( p, d, params, i, An, Bn, Cn, def, d.range=c(0,200), global, lwd=2, mgp=NULL, ... ){
# Plots the z function vs. d[i] in state i
d.seq <- seq(d.range[1], d.range[2],length.out=1000)
# The x-values
y.vals <- t( sapply( d.seq, function(d.i){
this.d <- d
this.d[i] <- d.i
return( zed_2(p, this.d, params, An, Bn, Cn, def )[i,] )
} ) )
# The y values
plot( d.seq, y.vals[,1], lwd=lwd, xlab='d', ylab='z', main=paste0( 'i = ', i ), type='l', mgp=mgp, ... )
# abline( 0, 1, lty=2 )
abline( v=d[i], lty=2 )
abline( h=p[i], lty=2 )
# The level
abline( v=global, lty=2, col='blue' )
# The global minimum
par(new = TRUE)
y.range <- c( max( min(c(0, y.vals[,2]), na.rm = TRUE), -10 ),
min( max(y.vals[,2], na.rm = TRUE ), 10 ) )
# The y limits for the derivative
plot( d.seq, y.vals[,2], type = "l", col='red', axes = FALSE, bty = "n",
xlab = "", ylab = "", ylim=y.range, mgp=mgp )
# Plot the derivative
axis(side=4, at = pretty(y.range), mgp=mgp-c(0,1,0) )
# mtext("z.p", side=4, line=3)
abline( h=1, lty=2, col='red' )
# Axis labelling
}
plot.q <- function( p, params, An=NULL, def=NULL ){
# Plots the q function for all the values of i
n <- length(params$R)
n.x <- ceiling( sqrt(n) )
n.y <- ceiling( n / n.x )
par(mfrow=c(n.y,n.x))
for( i in 1:n )
plot.q.i( p, params, i, An, def )
par(mfrow=c(1,1))
}
plot.q.i <- function( p, params, i, An=NULL, def=NULL ){
# Plots the q function vs. p[i] in state i
p.seq <- seq(0,1,by=.001)
# The x-values
n <- length(params$R)
# Number of states
if( is.null(An) ) An <- c(0)
if( is.null(def) ) def <- matrix(0,1,1)
# Fill in An and def if required
y.vals <- sapply( p.seq, function(p.i){
this.p <- p
this.p[i] <- p.i
return( q_fn(params$R, this.p, params$trans, params$lambda, params$phi, n,
params$cont.type, params$G, An, def ) )
} )
# The y values
plot( range(p.seq), range(y.vals), xlab='p', ylab='q', main=paste0( 'i = ', i ), type='n', ylim=c(0,1) )
for( j in 1:n ){
if( j == i ){
lines( p.seq, y.vals[j,], lwd=2, col='red' )
}else{
lines( p.seq, y.vals[j,], lwd=.5, col='black' )
}
}
# abline( h = ( 1 - params$lambda ) * params$phi / params$R[i], lty=2, lwd=1 )
abline( h = 1 / params$R[i], lty=2, lwd=1 )
abline( v = p[i], lty=2, lwd=1 )
legend( 'topright', paste0( 'R=', params$R[i]-1, '\nG=', params$G[i]-1), bty='n' )
}
plot.sol <- function( sol, ... ){
# Plots the outer solution
par(mfrow=c(2,2))
# Multiple plots
n <- length(params$R)
# Plot Q
plot( range(sol$d.grid), range(sol$Q), type='n', xlab='Debt value', ylab='',
main='Debt price', ... )
abline(v=sol$d.bar, lty=2, lwd=2, col=1:n )
abline(v=sol$d.grid, lty=1, lwd=.25 )
abline(h=c(0,1))
for( i in 1:n ) lines( sol$d.grid, sol$Q[i,], col=i, lwd=2 )
# Plot QE
plot( range(sol$d.grid), range(sol$QE), type='n', xlab='Debt value', ylab='',
main='Expected continuation debt price', ... )
abline(v=sol$d.bar, lty=2, lwd=2, col=1:n )
abline(v=sol$d.grid, lty=1, lwd=.25 )
abline(h=c(0,1))
for( i in 1:n ) lines( sol$d.grid, sol$QE[i,], col=i, lwd=2 )
# Plot P
plot( range(sol$d.grid), range(sol$P), type='n', xlab='Debt value', ylab='',
main='Default probability', ... )
abline(v=sol$d.bar, lty=2, lwd=2, col=1:n )
abline(v=sol$d.grid, lty=1, lwd=.25 )
abline(h=c(0,1))
for( i in 1:n ) lines( sol$d.grid, sol$P[i,], col=i, lwd=2 )
# Plot D.prime
plot( range(sol$d.grid), range(sol$D.prime), type='n', xlab='Debt value', ylab='',
main='Expected continuation debt', ... )
abline(v=sol$d.bar, lty=2, lwd=2, col=1:n )
abline(h=sol$d.bar, lty=2, col=1:n )
abline(v=sol$d.grid, lty=1, lwd=.25 )
abline(0,1,lty=2)
for( i in 1:n ) lines( sol$d.grid, sol$D.prime[i,], col=i, lwd=2 )
par(mfrow=c(1,1))
}
plot.err <- function(sol.err){
# Plots the errors from a solution
par(mfrow=c(2,2))
plot.err.i( sol.err$P, main='Default probability' )
plot.err.i( sol.err$Q, main='Debt price' )
plot.err.i( sol.err$D.prime, main='Ave. continuation debt' )
tb <- sapply( c('P','Q','D.prime'), function(x) c('Ave err'=mean(sol.err[[x]]),
'Ave abs err'=mean(abs(sol.err[[x]])) ) )
# Table of ave and ave abs errors
d.bar.err <- sol.err$d.bar.range - sol.err$d.bar
# The errors
nn <- length(params$R)
# Number of states
plot( c(1,nn), range(d.bar.err), type='n', xlab='State #', ylab='Upper-lower bound range',
xaxt='n', main='Inner vs. outer bounds check')
axis(1,at=1:nn,labels=format(1:nn,digits = 1))
abline(h=0,lwd=.5)
segments(1:nn,d.bar.err[,1],1:nn,d.bar.err[,2])
epsilon <- 0.02
segments(1:nn-epsilon,d.bar.err[,1],1:nn+epsilon,d.bar.err[,1])
segments(1:nn-epsilon,d.bar.err[,2],1:nn+epsilon,d.bar.err[,2])
# The consistency check for debt limits
par(mfrow=c(1,1))
return(tb)
}
plot.err.i <- function(err, ... ){
# Plots an inidividual error distribution
plot( density( abs(c(err)), from=0 ), xlab='Absolute error', ylab='Density', lwd=2, ... )
abline( v = quantile( abs(c(err)), c(.5, .75, .9, .95, .99)), lty=2 )
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.