diff --git a/src/QuantEcon.jl b/src/QuantEcon.jl index 85382ce1..ddbd3b54 100644 --- a/src/QuantEcon.jl +++ b/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, diff --git a/src/markov/markov_approx.jl b/src/markov/markov_approx.jl index 17b616a5..5f8b1d36 100644 --- a/src/markov/markov_approx.jl +++ b/src/markov/markov_approx.jl @@ -232,7 +232,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 +264,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 +321,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 +336,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 +358,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. (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 +419,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}, ::VAREstimationMethod) = nothing + + """ return standerdized AR(1) representation @@ -508,6 +532,22 @@ 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' """ construct one-dimensional evenly spaced grid of states @@ -526,7 +566,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 +580,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 +607,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 = qnwnorm(Nm, 0, 1) + y1D = repmat(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 +752,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]) """ diff --git a/test/test_markov_approx.jl b/test/test_markov_approx.jl index 1427ecb2..ad1168b8 100644 --- a/test/test_markov_approx.jl +++ b/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