Add quadrature method

REQUIRE

@@ -7,3 +7,4 @@ Compat 0.18.0
 StatsBase
 Optim 0.9.0
 NLopt 0.3.4
+FastGaussQuadrature 0.0.1

src/QuantEcon.jl

@@ -50,7 +50,7 @@ export
  simulate, simulate!, simulate_indices, simulate_indices!,
  period, is_irreducible, is_aperiodic, recurrent_classes,
  communication_classes, n_states,
- discrete_var, Even, Quantile,
+ discrete_var, Even, Quantile, Quadrature,
 
 # modeltools
  AbstractUtility, LogUtility, CRRAUtility, CFEUtility, EllipticalUtility, derivative,

src/markov/markov_approx.jl

@@ -13,6 +13,7 @@ https://lectures.quantecon.org/jl/finite_markov.html
 
 import Optim
 import NLopt
+import FastGaussQuadrature: gausshermite
 import Distributions: pdf, Normal, quantile
 
 std_norm_cdf(x::T) where {T <: Real} = 0.5 * erfc(-x/sqrt(2))
@@ -232,7 +233,8 @@ types specifying the method for `discrete_var`
 abstract type VAREstimationMethod end
 struct Even <: VAREstimationMethod end
 struct Quantile <: VAREstimationMethod end
-# immutable type Quadrature <: VAREstimationMethod end
+struct Quadrature <: VAREstimationMethod end
+
 
 doc"""
 
@@ -263,9 +265,8 @@ P, X = discrete_var(b, B, Psi, Nm, n_moments, method, n_sigmas)
 
 - `n_moments::Integer` : Desired number of moments to match. The default is 2.
 - `method::VAREstimationMethod` : Specify the method used to determine the grid
- points. Accepted inputs are `Even()` or `Quantile()`.
- Please see the paper for more details.
- NOTE: `Quadrature()` is not supported now.
+ points. Accepted inputs are `Even()`, `Quantile()`,
+ or `Quadrature()`. Please see the paper for more details.
 - `n_sigmas::Real` : If the `Even()` option is specified, `n_sigmas` is used to
  determine the number of unconditional standard deviations
  used to set the endpoints of the grid. The default is
@@ -321,6 +322,9 @@ function discrete_var(b::Union{Real, AbstractVector},
  error("n_moments must be either 1 or a positive even integer")
  end
 
+ # warning about persistency
+ warn_persistency(B, method)
+
  # Compute polynomial moments of standard normal distribution
  gaussian_moment = zeros(n_moments)
  c = 1
@@ -333,7 +337,7 @@ function discrete_var(b::Union{Real, AbstractVector},
  A, C, mu, Sigma = standardize_var(b, B, Psi, M)
 
  # Construct 1-D grids
- y1D, y1Dbounds = construct_1D_grid(Sigma, Nm, M, n_sigmas, method)
+ y1D, y1Dhelper = construct_1D_grid(Sigma, Nm, M, n_sigmas, method)
 
  # Construct all possible combinations of elements of the 1-D grids
  D = allcomb3(y1D')'
@@ -355,7 +359,7 @@ function discrete_var(b::Union{Real, AbstractVector},
  for ii = 1:(Nm^M)
 
  # Construct prior guesses for maximum entropy optimizations
- q = construct_prior_guess(cond_mean[:, ii], Nm, y1D, y1Dbounds, method)
+ q = construct_prior_guess(cond_mean[:, ii], Nm, y1D, y1Dhelper, method)
 
  # Make sure all elements of the prior are stricly positive
  q[q.<kappa] = kappa
@@ -371,17 +375,17 @@ function discrete_var(b::Union{Real, AbstractVector},
  # Maximum entropy optimization
  if n_moments == 1 # match only 1 moment
  temp[:, jj], _, _ = discrete_approximation(y1D[jj, :],
- X -> (reshape(X, 1, Nm)-cond_mean[jj, ii])/scaling_factor[jj],
+ X -> (X'-cond_mean[jj, ii])/scaling_factor[jj],
  [0.0], q[jj, :], [0.0])
  else # match 2 moments first
  p, lambda, moment_error = discrete_approximation(y1D[jj, :],
  X -> polynomial_moment(X, cond_mean[jj, ii], scaling_factor[jj], 2),
  [0; 1]./(scaling_factor[jj].^(1:2)), q[jj, :], lambda_guess)
- if norm(moment_error) > 1e-5 # if 2 moments fail, just match 1 moment
+ if !(norm(moment_error) < 1e-5) # if 2 moments fail, just match 1 moment
  warn("Failed to match first 2 moments. Just matching 1.")
  temp[:, jj], _, _ = discrete_approximation(y1D[jj, :],
-  X -> (X-cond_mean[jj,ii])/scaling_factor[jj],
- 0, q[jj, :], 0)
+ X -> (X'-cond_mean[jj, ii])/scaling_factor[jj],
+ [0.0], q[jj, :], [0.0])
  lambda_bar[(jj-1)*2+1:jj*2, ii] = zeros(2,1)
  elseif n_moments == 2
  lambda_bar[(jj-1)*2+1:jj*2, ii] = lambda
@@ -416,6 +420,27 @@ function discrete_var(b::Union{Real, AbstractVector},
  return MarkovChain(P, [X[:, i] for i in 1:Nm^M])
 end
 
+"""
+check persistency when `method` is `Quadrature` and give warning if needed
+
+##### Arguments
+- `B::Union{Real, AbstractMatrix}` : impact coefficient
+- `method::VAREstimationMethod` : method for grid making
+
+##### Returns
+- `nothing`
+
+"""
+function warn_persistency(B::AbstractMatrix, ::Quadrature)
+ if any(eig(B)[1] .> 0.9)
+ warn("The quadrature method may perform poorly for highly persistent processes.")
+ end
+ return nothing
+end
+warn_persistency(B::Real, method::Quadrature) = warn_persistency(fill(B, 1, 1), method)
+warn_persistency(::Union{Real, AbstractMatrix}, ::Union{Even, Quantile}) = nothing
+
+
 """
 
 return standerdized AR(1) representation
@@ -508,6 +533,23 @@ construct_prior_guess(cond_mean::AbstractVector, Nm::Integer,
  ::AbstractMatrix, y1Dbounds::AbstractMatrix, method::Quantile) =
  cdf.(Normal.(repmat(cond_mean, 1, Nm), 1), y1Dbounds[:, 2:end]) -
  cdf.(Normal.(repmat(cond_mean, 1, Nm), 1), y1Dbounds[:, 1:end-1])
+
+"""
+construct prior guess for quadrature grid method
+
+##### Arguments
+
+- `cond_mean::AbstractVector` : conditional Mean of each variable
+- `Nm::Integer` : number of grid points
+- `y1D::AbstractMatrix` : grid of variable
+- `weights::AbstractVector` : weights of grid `y1D`
+- `method::Quadrature` : method for grid making
+
+"""
+construct_prior_guess(cond_mean::AbstractVector, Nm::Integer,
+ y1D::AbstractMatrix, weights::AbstractVector, method::Quadrature) =
+ (pdf.(Normal.(repmat(cond_mean, 1, Nm), 1), y1D) ./ pdf.(Normal(0, 1), y1D)).*
+ (weights'/sqrt(pi))
 """
 
 construct one-dimensional evenly spaced grid of states
@@ -526,7 +568,7 @@ construct one-dimensional evenly spaced grid of states
 - `nothing` : `nothing` of type `Void`
 
 """
-function construct_1D_grid(Sigma::ScalarOrArray, Nm::Integer,
+function construct_1D_grid(Sigma::Union{Real, AbstractMatrix}, Nm::Integer,
  M::Integer, n_sigmas::Real, method::Even)
  min_sigmas = sqrt(minimum(eigfact(Sigma).values))
  y1Drow = collect(linspace(-min_sigmas*n_sigmas, min_sigmas*n_sigmas, Nm))'
@@ -540,7 +582,7 @@ construct one-dimensional quantile grid of states
 
 ##### Argument
 
-- `Sigma::ScalarOrArray` : variance-covariance matrix of the standardized process
+- `Sigma::AbstractMatrix` : variance-covariance matrix of the standardized process
 - `Nm::Integer` : number of grid points
 - `M::Integer` : number of variables (`M=1` corresponds to AR(1))
 - `n_sigmas::Real` : number of standard error determining end points of grid
@@ -567,6 +609,31 @@ construct_1D_grid(Sigma::Real, Nm::Integer,
 
 """
 
+construct one-dimensional quadrature grid of states
+
+##### Argument
+
+- `::ScalarOrArray` : not used
+- `Nm::Integer` : number of grid points
+- `M::Integer` : number of variables (`M=1` corresponds to AR(1))
+- `n_sigmas::Real` : not used
+- `method::Quadrature` : method for grid making
+
+##### Return
+
+- `y1D` : `M x Nm` matrix of variable grid
+- `weights` : weights on each grid
+
+"""
+function construct_1D_grid(::ScalarOrArray, Nm::Integer,
+ M::Integer, ::Real, method::Quadrature)
+ nodes, weights = gausshermite(Nm)
+ y1D = repmat(sqrt(2)*nodes', M, 1)
+ return y1D, weights
+end
+
+"""
+
 Return combinations of each column of matrix `A`.
 It is simiplifying `allcomb2` by using `gridmake` from QuantEcon
 
@@ -687,6 +754,9 @@ function discrete_approximation(D::AbstractVector, T::Function, TBar::AbstractVe
  moment_error = grad/minimum_value
  return p, lambda_bar, moment_error
 end
+# discrete_approximation(D::AbstractVector, T::Function, TBar::Real,
+# q::AbstractVector=ones(N)/N, lambda0::Real=0) =
+# discrete_approximation(D, T, [TBar], q, [lambda0])
 
 """
 

test/test_markov_approx.jl

@@ -71,8 +71,7 @@
  sigma2 = 1; # conditional variance
 
  mc = discrete_var(0, rho, sigma2, N, nMoments, Even())
- # `method` can be `Even()` or 'Quantile'
- # `Quadrature()` is not supported now
+ # `method` can be `Even()`, `Quantile()`, or `Quadrature()`
 
  # P1 computed by matlab code
  P1_matlab =[
@@ -117,14 +116,14 @@
  0.00121959272777097 4.11002290565958e-05 6.86068801631939e-06 1.84676508295551e-06 6.19295049673521e-07 2.31378890793673e-07 8.95788781706096e-08 0.0190492544153055 0.979680404921949]
 
  # D1 computed by matlab code
- D1_matlab = [-11.2940287556214 
- -6.85786965529225 
- -4.17854136278901 
- -2.00057722402480 
- 0 
- 2.00057722402480 
- 4.17854136278901 
- 6.85786965529225 
+ D1_matlab = [-11.2940287556214
+ -6.85786965529225
+ -4.17854136278901
+ -2.00057722402480
+ 0
+ 2.00057722402480
+ 4.17854136278901
+ 6.85786965529225
  11.2940287556214]
 
  # test transition matrix
@@ -382,9 +381,35 @@
  D2_matlab = [
  -0.0848530643859206 -0.0776879932266946 -0.0733604708271284 -0.0698427277337100 -0.0666114914670483 -0.0633802552003867 -0.0598625121069682 -0.0555349897074021 -0.0483699185481761 -0.0586888690883867 -0.0515237979291607 -0.0471962755295945 -0.0436785324361761 -0.0404472961695144 -0.0372160599028527 -0.0336983168094343 -0.0293707944098682 -0.0222057232506421 -0.0428863545423998 -0.0357212833831737 -0.0313937609836076 -0.0278760178901892 -0.0246447816235275 -0.0214135453568658 -0.0178958022634474 -0.0135682798638813 -0.00640320870465522 -0.0300408551387773 -0.0228757839795513 -0.0185482615799851 -0.0150305184865667 -0.0117992822199051 -0.00856804595324338 -0.00505030285982496 -0.000722780460258815 0.00644229069896724 -0.0182415729188723 -0.0110765017596462 -0.00674897936008010 -0.00323123626666167 0 0.00323123626666168 0.00674897936008010 0.0110765017596462 0.0182415729188723 -0.00644229069896724 0.000722780460258818 0.00505030285982496 0.00856804595324338 0.0117992822199051 0.0150305184865667 0.0185482615799852 0.0228757839795513 0.0300408551387773 0.00640320870465521 0.0135682798638813 0.0178958022634474 0.0214135453568658 0.0246447816235275 0.0278760178901892 0.0313937609836076 0.0357212833831737 0.0428863545423998 0.0222057232506421 0.0293707944098682 0.0336983168094343 0.0372160599028527 0.0404472961695144 0.0436785324361761 0.0471962755295945 0.0515237979291607 0.0586888690883867 0.0483699185481760 0.0555349897074021 0.0598625121069682 0.0633802552003867 0.0666114914670483 0.0698427277337100 0.0733604708271284 0.0776879932266946 0.0848530643859206
  -0.145665731585134 -0.0668724078570711 -0.0192832261459590 0.0194009214629551 0.0549343920375071 0.0904678626120590 0.129152010220973 0.176741191932085 0.255534515660148 -0.167243302213944 -0.0884499784858813 -0.0408607967747693 -0.00217664916585514 0.0333568214086968 0.0688902919832488 0.107574439592163 0.155163621303275 0.233956945031338 -0.180275611056287 -0.101482287328224 -0.0538931056171120 -0.0152089580081979 0.0203245125663541 0.0558579831409060 0.0945421307498201 0.142131312460932 0.220924636188995 -0.190869274170785 -0.112075950442722 -0.0644867687316099 -0.0258026211226958 0.00973084945185619 0.0452643200264082 0.0839484676353223 0.131537649346434 0.210330973074497 -0.200600123622641 -0.121806799894578 -0.0742176181834661 -0.0355334705745519 0 0.0355334705745520 0.0742176181834661 0.121806799894578 0.200600123622641 -0.210330973074497 -0.131537649346434 -0.0839484676353223 -0.0452643200264081 -0.00973084945185620 0.0258026211226958 0.0644867687316099 0.112075950442722 0.190869274170785 -0.220924636188995 -0.142131312460932 -0.0945421307498201 -0.0558579831409060 -0.0203245125663541 0.0152089580081979 0.0538931056171120 0.101482287328224 0.180275611056287 -0.233956945031338 -0.155163621303275 -0.107574439592163 -0.0688902919832488 -0.0333568214086968 0.00217664916585515 0.0408607967747693 0.0884499784858813 0.167243302213944 -0.255534515660148 -0.176741191932085 -0.129152010220973 -0.0904678626120590 -0.0549343920375071 -0.0194009214629551 0.0192832261459590 0.0668724078570711 0.145665731585134]
-
  @test isapprox(mc.p, P2_matlab, atol = 1e-6, rtol = 1e-6)
  # test state values
  @test all(isapprox.(mc.state_values, [D2_matlab[:, i] for i in 1:81], atol = 1e-6, rtol = 1e-6))
 
+ # test quadrature method
+ mc = discrete_var(0, rho, sigma2, N, nMoments, Quadrature())
+ P1_matlab = [
+ 0.965928500476864 0.0335561623394296 0.000511914985154753 3.41282320893676e-06 9.36593029189716e-09 9.27504215688250e-12 1.19716192091492e-13 1.57960519477171e-14 1.51230488749136e-15
+ 0.261803464438521 0.456272787943616 0.232464370220424 0.0456428985204484 0.00369588965275635 0.000119088911115234 1.32665822208226e-06 5.31037906465840e-08 1.20551105563997e-07
+ 0.0289335782187421 0.245371050107498 0.431199011778525 0.241697993985803 0.0491107678689957 0.00360255633129313 8.45663973919190e-05 4.65221553851222e-07 1.00901972676700e-08
+ 0.00129342416375824 0.0429398178267666 0.244970263155592 0.412007463325491 0.243253715549902 0.0518252887496408 0.00364496551253315 6.49233683356058e-05 1.38347980459142e-07
+ 2.23458474109999e-05 0.00278914152743075 0.0499164081274372 0.244097503218571 0.406349203242795 0.244097502746911 0.0499164079316752 0.00278914151054822 2.23458472205773e-05
+ 1.38347980459062e-07 6.49233683355849e-05 0.00364496551253256 0.0518252887496380 0.243253715549904 0.412007463325499 0.244970263155590 0.0429398178267626 0.00129342416375791
+ 1.00901961081846e-08 4.65221519912425e-07 8.45663936588306e-05 0.00360255624680754 0.0491107674227084 0.241697993895957 0.431199012942789 0.245371049947396 0.0289335778389678
+ 1.20551105544828e-07 5.31037906407268e-08 1.32665822198243e-06 0.000119088911109472 0.00369588965265401 0.0456428985198852 0.232464370219936 0.456272787944799 0.261803464438499
+ 1.51230505234139e-15 1.57960534189186e-14 1.19716201601543e-13 9.27504277505984e-12 9.36593079980241e-09 3.41282335162177e-06 0.000511915000010395 0.0335561628535355 0.965928499947760]
+ D1_matlab = [
+ -4.51274586339979
+ -3.20542900285647
+ -2.07684797867783
+ -1.02325566378913
+ -5.69133976269097e-17
+ 1.02325566378913
+ 2.07684797867783
+ 3.20542900285647
+ 4.51274586339979]
+
+ # test transition matrix
+ @test isapprox(mc.p, P1_matlab, atol = 1e-7, rtol = 1e-7)
+ # test state values
+ @test isapprox(mc.state_values, D1_matlab)
 end # @testset