context('Checking NGM model')
test_that("Integration", {
params <- list( A=1, alpha=1/3, delta=.05, gamma=1, delta=.02, betta=.95 )
upper <- c( 1, 1.5 )
lower <- c( -.5, 0 )
coeff <- matrix( c( 1, .1, .15 ), 3, 1)
exog <- matrix( 0, 1, 1 )
endog <- matrix( 1, 2, 1 )
# Variables used
rho <- .8
sig.eps <- .1
n.mc <- 100000 # 1000000
# Number of Monte Carlo shocks. Works better & slower with bigger numbers.
n.quad <- 5
# Number of quadrature nodes
set.seed(12345)
exog.innov <- matrix( rnorm( n.mc, 0, sig.eps ), n.mc, 1 )
# Monte Carlo integration shocks
err.mc <- euler_hat_ngm( exog, endog, exog.innov, params, coeff, 1, 1, rho,
n.mc, 1, upper, lower, FALSE, rep(1,n.mc)/n.mc )
nodes <- quad_nodes_1d( n.quad, sig.eps, 0 )
wts <- quad_weights_1d( n.quad )
err.quad <- euler_hat_ngm( exog, endog, nodes, params,
coeff, 1, 1, rho, n.quad, 1,
upper, lower, FALSE, wts )
expect_true( abs( err.mc - err.quad ) < 1e-04 )
# Check that the integrals are close(ish). MCintegration is quite bad!
})
test_that("Derivatives", {
skip("Derivatives not functional at present")
params <- list( A=1, alpha=1/3, delta=.05, gamma=1, delta=.02, betta=.95 )
upper <- c( 1, 1.5 )
lower <- c( -.5, 0 )
N <- 2
# approximation order
coeff <- matrix( c( 1, .1, .15, .1, .1, -.05 ), 6, 1)
exog <- 0
endog <- 1
exog_shk <- .01
# Variables used
rho <- .8
sig.eps <- .1
n.quad <- 7
# Number of quadrature nodes
##### The integrand #####
## Ordinary polynomials
base <- integrand_ngm( exog, endog, exog_shk, params, coeff, 1, 1, 2, upper,
lower, FALSE)
deriv.x <- integrand_ngm_D( exog, endog, exog_shk, params, coeff, 1, 1, 2,
upper, lower, FALSE)
# The exact derivative
deriv.fd <- 0 * deriv.x
# Initialize the first-diff approx
inc <- 1e-06 * diag( 6 )
# The matrix of increments
for( i in 1:6 )
deriv.fd[i] <- 1e06 * ( integrand_ngm( exog, endog, exog_shk, params,
coeff + inc[i,],
1, 1, 2, upper, lower, FALSE) - base )
# Fill in the vector of first derivatives
expect_equal( deriv.x, deriv.fd, tolerance=1e-04 )
# Need to have low tolerance because FD is only an approximation
## Chebychev polynomials
base <- integrand_ngm( exog, endog, exog_shk, params, coeff, 1, 1, 2, upper,
lower, TRUE )
deriv.x <- integrand_ngm_D( exog, endog, exog_shk, params, coeff, 1, 1, 2,
upper, lower, TRUE )
# The exact derivative
deriv.fd <- 0 * deriv.x
# Initialize the first-diff approx
inc <- 1e-06 * diag( 6 )
# The matrix of increments
for( i in 1:6 )
deriv.fd[i] <- 1e06 * ( integrand_ngm( exog, endog, exog_shk, params,
coeff + inc[i,],
1, 1, 2, upper, lower, TRUE) - base )
# Fill in the vector of first derivatives
expect_equal( deriv.x, deriv.fd, tolerance=1e-05 )
# Need to have low tolerance because FD is only an approximation
##### The integral #####
## Ordinary polynomials
exog <- matrix(0,1,1)
nodes <- quad_nodes_1d( n.quad, sig.eps, 0 )
wts <- quad_weights_1d( n.quad )
base <- err_ngm( exog, endog, nodes, params, coeff, 1, 1, rho,
n.quad, 2, upper, lower, FALSE, wts )
# Case where the sign of the error derivative is important
deriv.x <- err_ngm_D( exog, endog, nodes, params, coeff, 1, 1, rho,
n.quad, 2, upper, lower, FALSE, wts )
# The exact derivative
deriv.fd <- 0 * deriv.x
# Initialize the first-diff approx
inc <- 1e-06 * diag( 6 )
# The matrix of increments
for( i in 1:6 )
deriv.fd[i] <- 1e06 * ( err_ngm( exog, endog, nodes, params,
coeff + inc[i,], 1, 1, rho,
n.quad, 2, upper, lower, FALSE, wts ) - base )
# Fill in the vector of first derivatives
expect_equal( deriv.x, deriv.fd, tolerance=1e-04 )
# Need to have low tolerance because FD is only an approximation
exog <- matrix(0.2,1,1)
base <- err_ngm( exog, endog, nodes, params, coeff, 1, 1, rho,
n.quad, 2, upper, lower, FALSE, wts )
# Case where the sign of the error derivative is not important
deriv.x <- err_ngm_D( exog, endog, nodes, params, coeff, 1, 1, rho,
n.quad, 2, upper, lower, FALSE, wts )
# The exact derivative
deriv.fd <- 0 * deriv.x
# Initialize the first-diff approx
inc <- 1e-06 * diag( 6 )
# The matrix of increments
for( i in 1:6 )
deriv.fd[i] <- 1e06 * ( err_ngm( exog, endog, nodes, params,
coeff + inc[i,], 1, 1, rho,
n.quad, 2, upper, lower, FALSE, wts ) - base )
# Fill in the vector of first derivatives
expect_equal( deriv.x, deriv.fd, tolerance=1e-04 )
# Need to have low tolerance because FD is only an approximation
## Chebychev polynomials
exog <- matrix(0,1,1)
nodes <- quad_nodes_1d( n.quad, sig.eps, 0 )
wts <- quad_weights_1d( n.quad )
base <- err_ngm( exog, endog, nodes, params, coeff, 1, 1, rho,
n.quad, 2, upper, lower, TRUE, wts )
# Case where the sign of the error derivative is important
deriv.x <- err_ngm_D( exog, endog, nodes, params, coeff, 1, 1, rho,
n.quad, 2, upper, lower, TRUE, wts )
# The exact derivative
deriv.fd <- 0 * deriv.x
# Initialize the first-diff approx
inc <- 1e-06 * diag( 6 )
# The matrix of increments
for( i in 1:6 )
deriv.fd[i] <- 1e06 * ( err_ngm( exog, endog, nodes, params,
coeff + inc[i,], 1, 1, rho,
n.quad, 2, upper, lower, TRUE, wts ) - base )
# Fill in the vector of first derivatives
expect_equal( deriv.x, deriv.fd, tolerance=1e-05 )
# Need to have low tolerance because FD is only an approximation
exog <- matrix(0.6,1,1)
base <- err_ngm( exog, endog, nodes, params, coeff, 1, 1, rho,
n.quad, 2, upper, lower, TRUE, wts )
# Case where the sign of the error derivative is not important
deriv.x <- err_ngm_D( exog, endog, nodes, params, coeff, 1, 1, rho,
n.quad, 2, upper, lower, TRUE, wts )
# The exact derivative
deriv.fd <- 0 * deriv.x
# Initialize the first-diff approx
inc <- 1e-06 * diag( 6 )
# The matrix of increments
for( i in 1:6 )
deriv.fd[i] <- 1e06 * ( err_ngm( exog, endog, nodes, params,
coeff + inc[i,], 1, 1, rho,
n.quad, 2, upper, lower, TRUE, wts ) - base )
# Fill in the vector of first derivatives
expect_equal( deriv.x, deriv.fd, tolerance=1e-04 )
# Need to have low tolerance because FD is only an approximation
} )
test_that("Checking full-depreciation solution", {
params.full.dep <- list( A=1, alpha=.3, delta=1, gamma=1, betta=.99,
rho=.0, sig.eps=.01 )
k.ss <- ( ( 1 / params.full.dep$betta - 1 + params.full.dep$delta ) /
( params.full.dep$A * params.full.dep$alpha ) ) ^ ( 1 / ( params.full.dep$alpha - 1 ) )
# The steady state
sd.x<- sqrt( params.full.dep$sig.eps / ( 1 - params.full.dep$rho ^ 2 ) )
upper <- c( 5 * sd.x, k.ss * 1.4 )
lower <- c( -5* sd.x, k.ss * 0.6 )
# Create the bounds
opt <- list( model='ngm', lags=1, n.exog=1, n.endog=1, N=1, cheby=FALSE,
upper = upper, lower=lower, quad=TRUE, n.quad=7,
err.tol=1e-05, n.iter=10, burn=0, kappa=1, n.sim=10000,
eps = .3, delta=.02, endog.init=k.ss )
# Solution options
a.b.range <- upper - lower
coeff.init <- matrix( c( k.ss,
params.full.dep$alpha * a.b.range[2] / 2,
k.ss * a.b.range[1] / 2 ), 3, 1)
# The linear-approx solution evaluated @ the steady state:
# k' = alpha * betta * exp(x) * A * k ^ alpha
# ~= k.ss + k.ss * x + alpha * ( k - k.ss )
# NB: The range is already centered on k.ss, so can just use
# these coeffs
k.ss.prime <- endog_update( exog = c(0), k.ss, coeff.init, 1, 1, 1, upper, lower, FALSE )
expect_true( k.ss == k.ss.prime )
k.ss.euler <-
integrand_ngm( matrix( 0, 1, 1), matrix( k.ss, nrow=2, ncol=1), 0,
params.full.dep, coeff.init, 1, 1, 1, upper, lower, FALSE )
expect_equal( k.ss.euler * params.full.dep$betta, 1 )
nodes <- quad_nodes_1d( opt$n.quad, params.full.dep$sig.eps, 0 )
wts <- quad_weights_1d( opt$n.quad )
single.integral <-
euler_hat_ngm( matrix( 0, 1, 1), matrix( k.ss, nrow=2, ncol=1), nodes,
params.full.dep, coeff.init, opt$n.exog, opt$n.endog,
params.full.dep$rho, opt$n.quad, opt$N, upper, lower, FALSE, wts )
expect_equal( k.ss, c( single.integral ), tolerance=4e-05 )
# These are not exactly equal because we integrate over a linear
# approximation to the solution.
mult.pts <-
euler_hat( coeff.init, matrix( c( 0, k.ss, 0, k.ss ), nrow=1 ), 'ngm', opt$lags,
params.full.dep, opt$n.exog, opt$n.endog, params.full.dep$rho,
params.full.dep$sig.eps, 0, opt$N, upper, lower, FALSE,
matrix(0,1,1,), TRUE, opt$n.quad )
expect_equal( single.integral, mult.pts )
exog.sim <- ar1_sim( opt$n.sim * opt$kappa + opt$burn,
params.full.dep$rho, params.full.dep$sig.eps )
endog.sim <- endog_sim( opt$n.sim, exog.sim, coeff.init, opt$N, opt$upper,
opt$lower, k.ss, FALSE, opt$kappa, opt$burn, opt$lags )
indices <- p_eps_cheap_const_idx( endog.sim[,3:4], opt$eps, opt$delta )
sum(indices)
X <- endog.sim[indices==1,]
# State reduction
should.be.coeff.init <-
coeff_reg( endog.sim[,2], endog.sim[,c(1,4)], opt$N,
opt$lower, opt$upper, opt$cheby )
expect_equal( coeff.init, should.be.coeff.init )
# Check that the regression does the right thing
coeff.new <- coeff.init
for( i in 1:5 ){
endog.sim <- endog_sim( opt$n.sim, exog.sim, coeff.new, opt$N, opt$upper,
opt$lower, k.ss, FALSE, opt$kappa, opt$burn, opt$lags )
# indices <- p_eps_cheap_const_idx( endog.sim[,3:4], opt$eps, opt$delta )
# sum(indices)
# X <- endog.sim[indices==1,]
k.hat <- euler_hat( coeff.new, endog.sim, 'ngm', opt$lags,
params.full.dep, opt$n.exog, opt$n.endog, params.full.dep$rho,
params.full.dep$sig.eps, 0, opt$N, upper, lower, FALSE,
matrix(0,1,1,), TRUE, opt$n.quad )
coeff.old <- coeff.new
coeff.new <- coeff_reg( k.hat, endog.sim[,c(1,4)], opt$N, opt$lower, opt$upper, opt$cheby )
}
## FIGURE OUT WHAT IS GOING ON HERE ##
# sol <- sol.iterate( coeff.init, opt, params.full.dep, TRUE )
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.