diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m index 66cec2f3f..7ff62baf6 100644 --- a/matlab/default_option_values.m +++ b/matlab/default_option_values.m @@ -256,7 +256,7 @@ particle.number_of_particles = 5000; % Set the default approximation order (Smolyak) particle.smolyak_accuracy = 3; % By default we don't use pruning -particle.pruning = 0; +particle.pruning = false; % Set the Gaussian approximation method for particles distributions particle.approximation_method = 'unscented'; % Set unscented parameters alpha, beta and kappa for gaussian approximation diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m index f3c70fc00..8455944a7 100644 --- a/matlab/non_linear_dsge_likelihood.m +++ b/matlab/non_linear_dsge_likelihood.m @@ -138,6 +138,7 @@ else ReducedForm.ghxx = dr.ghxx(restrict_variables_idx,:); ReducedForm.ghuu = dr.ghuu(restrict_variables_idx,:); ReducedForm.ghxu = dr.ghxu(restrict_variables_idx,:); + ReducedForm.ghs2 = dr.ghs2(restrict_variables_idx,:); if DynareOptions.order==3 ReducedForm.ghxxx = dr.ghxxx(restrict_variables_idx,:); ReducedForm.ghuuu = dr.ghuuu(restrict_variables_idx,:); diff --git a/matlab/particles b/matlab/particles index bb16df7f2..4a679b26e 160000 --- a/matlab/particles +++ b/matlab/particles @@ -1 +1 @@ -Subproject commit bb16df7f2e6ddbc81e23b9661ee100f5e8dc675a +Subproject commit 4a679b26eef0463169705cde460079feddc14819 diff --git a/mex/sources/local_state_space_iterations/local_state_space_iteration_3.f08 b/mex/sources/local_state_space_iterations/local_state_space_iteration_3.f08 index 588b495c2..f235fae49 100644 --- a/mex/sources/local_state_space_iterations/local_state_space_iteration_3.f08 +++ b/mex/sources/local_state_space_iterations/local_state_space_iteration_3.f08 @@ -24,13 +24,15 @@ module pparticle_3 type tdata_3 integer :: n, m, s, q, numthreads, xx_size, uu_size, xxx_size, uuu_size - real(real64), pointer, contiguous :: yhat3(:,:), & - &e(:,:), ghx(:,:), ghu(:,:), constant(:), ghxu(:,:), ghxx(:,:), & - &ghuu(:,:), ghxxx(:,:), ghuuu(:,:), ghxxu(:,:), ghxuu(:,:), ghxss(:,:), & - &ghuss(:,:), ss(:), y3(:,:) - real(real64), pointer :: yhat2(:,:), yhat1(:,:), y2(:,:), y1(:,:) + real(real64), pointer, contiguous :: e(:,:), ghx(:,:), ghu(:,:), & + &ghxu(:,:), ghxx(:,:), ghuu(:,:), ghs2(:), & + &ghxxx(:,:), ghuuu(:,:), ghxxu(:,:), ghxuu(:,:), ghxss(:,:), ghuss(:,:), & + &ss(:), y3(:,:) + real(real64), pointer :: yhat3(:,:), yhat2(:,:), yhat1(:,:), ylat3(:,:), & + &ylat2(:,:), ylat1(:,:) type(index), pointer, contiguous :: xx_idcs(:), uu_idcs(:), & &xxx_idcs(:), uuu_idcs(:) + integer, pointer, contiguous :: xx_nbeq(:) end type tdata_3 type(tdata_3) :: td3 @@ -66,7 +68,7 @@ contains do is=start,end do im=1,td3%m ! y3 = ybar + ½ghss - td3%y3(im,is) = td3%constant(im) + td3%y3(im,is) = td3%ss(im)+0.5_real64*td3%ghs2(im) ! y3 += ghx·ŷ+(3/6)·ghxss·ŷ + first n folded indices for ½ghxx·ŷ⊗ŷ ! + first n folded indices for (1/6)ghxxx·ŷ⊗ŷ⊗ŷ do j=1,td3%n @@ -153,19 +155,18 @@ contains end subroutine thread_eval_3 - ! Fills y1 and y2 as - ! y1 = ybar + ghx·ŷ1 + ghu·ε - ! y2 = ybar + ½ghss + ghx·ŷ2 + ghu·ε + ½ghxx·ŷ1⊗ŷ1 + ½ghuu·ε⊗ε + ghxu·ŷ1⊗ε - ! y3 = ybar + ghx·ŷ3 + ghu·ε + ghxx·ŷ1⊗ŷ2 + ghuu·ε⊗ε + ghxu·ŷ1⊗ε + ghxu·ŷ2⊗ε - ! + (1/6)·ghxxx·ŷ1⊗ŷ1⊗ŷ1 + (1/6)·ghuuu·ε⊗ε⊗ε + (3/6)·ghxxu·ŷ1⊗ŷ1⊗ε - ! + (3/6)·ghxuu·ŷ1⊗ε⊗ε + (3/6)·ghxss·ŷ1 + (3/6)·ghuss·ε + ! Fills ylat1, ylat2, ylat3 and y3 as + ! ylat1 = ghx·ŷ1 + ghu·ε + ! ylat2 = ½ghss + ghx·ŷ2 + ½ghxx·ŷ1⊗ŷ1 + ½ghuu·ε⊗ε + ghxu·ŷ1⊗ε + ! ylat3 = ghx·ŷ3 + ghxx·ŷ1⊗ŷ2 + ghxu·ŷ2⊗ε + (1/6)·ghxxx·ŷ1⊗ŷ1⊗ŷ1 + ! + (1/6)·ghuuu·ε⊗ε⊗ε + (3/6)·ghxxu·ŷ1⊗ŷ1⊗ε + ! + (3/6)·ghxuu·ŷ1⊗ε⊗ε + (3/6)·ghxss·ŷ1 + (3/6)·ghuss·ε + ! y3 = ybar + ylat1 + ylat2 + ylat3 ! in td3 subroutine thread_eval_3_pruning(arg) bind(c) type(c_ptr), intent(in), value :: arg integer, pointer :: ithread - integer :: is, im, j, k, start, end, q, r - real(real64) :: x, y - + integer :: is, im, j, k, start, end, q, r, j1, j2 ! Checking that the thread number got passed as argument if (.not. c_associated(arg)) then call mexErrMsgTxt("No argument was passed to thread_eval") @@ -184,79 +185,75 @@ contains do is=start,end do im=1,td3%m - ! y1 = ybar - ! y2 = ybar + ½ghss - ! y3 = ybar - td3%y1(im,is) = td3%ss(im) - td3%y2(im,is) = td3%constant(im) - td3%y3(im,is) = td3%ss(im) + ! y1 = 0 + ! y2 = ½ghss + ! y3 = 0 + td3%ylat1(im,is) = td3%ss(im) + td3%ylat2(im,is) = td3%ss(im)+0.5_real64*td3%ghs2(im) + td3%ylat3(im,is) = td3%ss(im) ! y1 += ghx·ŷ1 ! y2 += ghx·ŷ2 + first n folded indices for ½ghxx·ŷ1⊗ŷ1 - ! y3 += ghx·ŷ3 +(3/6)·ghxss·ŷ - ! + first n folded indices for ghxx·ŷ1⊗ŷ2 + ! y3 += ghx·ŷ3 +(3/6)·ghxss·ŷ1 ! + first n folded indices for (1/6)ghxxx·ŷ1⊗ŷ1⊗ŷ1 do j=1,td3%n - td3%y1(im,is) = td3%y1(im,is) + td3%ghx(j,im)*td3%yhat1(j,is) - td3%y2(im,is) = td3%y2(im,is) + td3%ghx(j,im)*td3%yhat2(j,is) +& + td3%ylat1(im,is) = td3%ylat1(im,is)+& + &td3%ghx(j,im)*td3%yhat1(j,is) + td3%ylat2(im,is) = td3%ylat2(im,is)+& + &td3%ghx(j,im)*td3%yhat2(j,is)+& &0.5_real64*td3%ghxx(j,im)*td3%yhat1(1, is)*td3%yhat1(j, is) - td3%y3(im,is) = td3%y3(im,is) + td3%ghx(j,im)*td3%yhat3(j,is) +& - &0.5_real64*td3%ghxss(j,im)*td3%yhat1(j,is) +& - &td3%ghxx(j,im)*td3%yhat1(1, is)*td3%yhat2(j, is)+& - &(1._real64/6._real64)*td3%ghxxx(j,im)*td3%yhat1(1, is)*& - &td3%yhat1(1, is)*td3%yhat1(j,is) + td3%ylat3(im,is) = td3%ylat3(im,is)+& + &td3%ghx(j,im)*td3%yhat3(j,is)+& + &0.5_real64*td3%ghxss(j,im)*td3%yhat1(j,is)+& + &(1._real64/6._real64)*td3%ghxxx(j,im)*& + &td3%yhat1(1, is)*td3%yhat1(1, is)*td3%yhat1(j,is) ! y2 += + ghxu·ŷ1⊗ε - ! y3 += + ghxu·ŷ1⊗ε + ghxu·ŷ2⊗ε + ! y3 += + ghxu·ŷ2⊗ε ! + first n*q folded indices of (3/6)·ghxxu·ŷ1⊗ŷ1⊗ε do k=1,td3%q - td3%y2(im,is) = td3%y2(im,is)+& + td3%ylat2(im,is) = td3%ylat2(im,is)+& &td3%ghxu(td3%q*(j-1)+k,im)*& &td3%yhat1(j, is)*td3%e(k, is) - td3%y3(im,is) = td3%y3(im,is)+& + td3%ylat3(im,is) = td3%ylat3(im,is)+& &td3%ghxu(td3%q*(j-1)+k,im)*& - &(td3%yhat1(j, is)+td3%yhat2(j, is))*td3%e(k, is)+& - &0.5_real64*td3%ghxxu(td3%q*(j-1)+k,im)*td3%yhat1(1, is)*& - &td3%yhat1(j, is)*td3%e(k, is) + &td3%yhat2(j, is)*td3%e(k, is)+& + &0.5_real64*td3%ghxxu(td3%q*(j-1)+k,im)*& + &td3%yhat1(1, is)*td3%yhat1(j, is)*td3%e(k, is) end do end do - ! y1 += +ghu·ε - ! y2 += +ghu·ε + first q folded indices for ½ghuu·ε⊗ε - ! y3 += +ghu·ε + first q folded indices for ghuu·ε⊗ε + ! y1 += + ghu·ε + ! y2 += + first q folded indices for ½ghuu·ε⊗ε + ! y3 += + (3/6)·ghuss·ε ! + first n*q folded indices of (3/6)·ghxuu·ŷ1⊗ε⊗ε ! + first n folded indices of (1/6)·ghuuu·ε⊗ε⊗ε do j=1,td3%q - x = td3%ghu(j,im)*td3%e(j,is) - y = td3%ghuu(j,im)*td3%e(1, is)*td3%e(j, is) - td3%y1(im,is) = td3%y1(im,is) + x - td3%y2(im,is) = td3%y2(im,is) + x + 0.5_real64*y - td3%y3(im,is) = td3%y3(im,is) + x + y +& - &td3%ghuss(j,im)*td3%e(j,is)+& - &(1._real64/6._real64)*td3%ghuuu(j,im)*td3%e(1, is)*td3%e(1, is)*& - &td3%e(j, is) + td3%ylat1(im,is) = td3%ylat1(im,is) + td3%ghu(j,im)*td3%e(j,is) + td3%ylat2(im,is) = td3%ylat2(im,is) + 0.5_real64*td3%ghuu(j,im)*td3%e(1, is)*td3%e(j, is) + td3%ylat3(im,is) = td3%ylat3(im,is)+& + &0.5_real64*td3%ghuss(j,im)*td3%e(j,is)+& + &(1._real64/6._real64)*td3%ghuuu(j,im)*& + &td3%e(1, is)*td3%e(1, is)*td3%e(j, is) do k=1,td3%n - td3%y3(im,is) = td3%y3(im,is) + & + td3%ylat3(im,is) = td3%ylat3(im,is)+& &0.5_real64*td3%ghxuu(td3%uu_size*(k-1)+j,im)*& &td3%yhat1(k, is)*td3%e(1, is)*td3%e(j, is) end do end do ! y2 += remaining ½ghxx·ŷ1⊗ŷ1 terms - ! y3 += remaining ghxx·ŷ1⊗ŷ2 terms - ! + the next terms starting from n+1 up to xx_size + ! y3 += + the next terms starting from n+1 up to xx_size ! of (1/6)ghxxx·ŷ1⊗ŷ1⊗ŷ1 - ! + remaining terms of (3/6)·ghxxu·ŷ1⊗ŷ1⊗ε + ! + remaining terms of (3/6)·ghxxu·ŷ1⊗ŷ1⊗ε do j=td3%n+1,td3%xx_size - td3%y2(im,is) = td3%y2(im,is) + & + td3%ylat2(im,is) = td3%ylat2(im,is)+& &0.5_real64*td3%ghxx(j,im)*& &td3%yhat1(td3%xx_idcs(j)%coor(1), is)*& &td3%yhat1(td3%xx_idcs(j)%coor(2), is) - td3%y3(im,is) = td3%y3(im,is) + & - &td3%ghxx(j,im)*& - &td3%yhat1(td3%xx_idcs(j)%coor(1), is)*& - &td3%yhat2(td3%xx_idcs(j)%coor(2), is)+& - &(1._real64/6._real64)*td3%ghxxx(j,im)*td3%yhat1(1, is)*& + td3%ylat3(im,is) = td3%ylat3(im,is)+& + &(1._real64/6._real64)*td3%ghxxx(j,im)*& + &td3%yhat1(1, is)*& &td3%yhat1(td3%xx_idcs(j)%coor(1), is)*& &td3%yhat1(td3%xx_idcs(j)%coor(2), is) - do k=1,td3%n - td3%y3(im,is) = td3%y3(im,is) + & + do k=1,td3%q + td3%ylat3(im,is) = td3%ylat3(im,is)+& &0.5_real64*td3%ghxxu(td3%q*(j-1)+k,im)*& &td3%yhat1(td3%xx_idcs(j)%coor(1), is)*& &td3%yhat1(td3%xx_idcs(j)%coor(2), is)*& @@ -264,21 +261,21 @@ contains end do end do ! y2 += remaining ½ghuu·ε⊗ε terms - ! y3 += remaining ghuu·ε⊗ε terms - ! + remaining terms of (3/6)·ghxuu·ŷ⊗ε⊗ε - ! + the next uu_size terms starting from q+1 - ! of (1/6)·ghuuu·ε⊗ε⊗ε + ! y3 += + remaining terms of (3/6)·ghxuu·ŷ⊗ε⊗ε + ! + the next uu_size terms starting from q+1 + ! of (1/6)·ghuuu·ε⊗ε⊗ε do j=td3%q+1,td3%uu_size - x = td3%ghuu(j,im)*& + td3%ylat2(im,is) = td3%ylat2(im,is)+& + &0.5_real64*td3%ghuu(j,im)*& &td3%e(td3%uu_idcs(j)%coor(1), is)*& &td3%e(td3%uu_idcs(j)%coor(2), is) - td3%y2(im,is) = td3%y2(im,is)+0.5_real64*x - td3%y3(im,is) = td3%y3(im,is)+x+& - &(1._real64/6._real64)*td3%ghuuu(j,im)*td3%e(1, is)*& + td3%ylat3(im,is) = td3%ylat3(im,is)+& + &(1._real64/6._real64)*td3%ghuuu(j,im)*& + &td3%e(1, is)*& &td3%e(td3%uu_idcs(j)%coor(1), is)*& &td3%e(td3%uu_idcs(j)%coor(2), is) do k=1,td3%n - td3%y3(im,is) = td3%y3(im,is) + & + td3%ylat3(im,is) = td3%ylat3(im,is)+& &0.5_real64*td3%ghxuu(td3%uu_size*(k-1)+j,im)*& &td3%yhat1(k, is)*& &td3%e(td3%uu_idcs(j)%coor(1), is)*& @@ -287,7 +284,7 @@ contains end do ! y3 += remaining (1/6)·ghxxx·ŷ⊗ŷ⊗ŷ terms do j=td3%xx_size+1,td3%xxx_size - td3%y3(im,is) = td3%y3(im,is)+& + td3%ylat3(im,is) = td3%ylat3(im,is)+& &(1._real64/6._real64)*td3%ghxxx(j,im)*& &td3%yhat1(td3%xxx_idcs(j)%coor(1), is)*& &td3%yhat1(td3%xxx_idcs(j)%coor(2), is)*& @@ -295,12 +292,32 @@ contains end do ! y3 += remaining (1/6)ghuuu·ε⊗ε⊗ε terms do j=td3%uu_size+1,td3%uuu_size - td3%y3(im,is) = td3%y3(im,is) + & + td3%ylat3(im,is) = td3%ylat3(im,is)+& &(1._real64/6._real64)*td3%ghuuu(j,im)*& &td3%e(td3%uuu_idcs(j)%coor(1), is)*& &td3%e(td3%uuu_idcs(j)%coor(2), is)*& &td3%e(td3%uuu_idcs(j)%coor(3), is) end do + ! y3 += first n folded indices for ghxx·ŷ1⊗ŷ2 + do j=1,td3%xx_size + j1 = td3%xx_idcs(j)%coor(1) + j2 = td3%xx_idcs(j)%coor(2) + if (j1 == j2) then + td3%ylat3(im,is) = td3%ylat3(im,is)+& + &td3%ghxx(j,im)*& + &td3%yhat1(j1,is)*& + &td3%yhat2(j2,is) + else + td3%ylat3(im,is) = td3%ylat3(im,is)+& + &0.5_real64*td3%ghxx(j,im)*(& + &td3%yhat1(j1,is)*& + &td3%yhat2(j2,is)+& + &td3%yhat1(j2,is)*& + &td3%yhat2(j1,is)) + end if + end do + td3%y3(im,is) = td3%ylat1(im,is)+td3%ylat2(im,is)+& + &td3%ylat3(im,is)-2*td3%ss(im) end do end do @@ -310,29 +327,34 @@ end module pparticle_3 ! The code of the local_state_space_iteration_3 routine ! Input: -! prhs[1] yhat3 [double] n×s array, time t particles. -! prhs[2] e [double] q×s array, time t innovations. -! prhs[3] ghx [double] m×n array, first order reduced form. -! prhs[4] ghu [double] m×q array, first order reduced form. -! prhs[5] constant [double] m×1 array, deterministic steady state + -! third order correction for the union of -! the states and observed variables. -! prhs[6] ghxx [double] m×n² array, second order reduced form. -! prhs[7] ghuu [double] m×q² array, second order reduced form. -! prhs[8] ghxu [double] m×nq array, second order reduced form. -! prhs[9] ghxxx [double] m×n array, third order reduced form. -! prhs[10] ghuuu [double] m×q array, third order reduced form. -! prhs[11] ghxxu [double] m×n²q array, third order reduced form. -! prhs[12] ghxuu [double] m×nq² array, third order reduced form. -! prhs[13] ghxss [double] m×n array, third order reduced form. -! prhs[14] ghuss [double] m×q array, third order reduced form. -! prhs[15] yhat2 [double] [OPTIONAL] 2n×s array, time t particles up to 2nd order pruning additional latent variables. The first n rows concern the pruning first-order latent variables, while the last n rows concern the pruning 2nd-order latent variables -! prhs[16] ss [double] [OPTIONAL] m×1 array, steady state for the union of the states and the observed variables (needed for the pruning mode). -! prhs[15 or 17] [double] num of threads +! prhs[1] yhat [double] n×s array, time t particles if pruning is false +! 3n×s array, time t particles if pruning is true +! Rows 1 to n contain the pruned first order. +! Rows n+1 to 2*n contain the pruned second order. +! Rows 2*n+1 to 3*n contain the pruned third order. +! prhs[2] e [double] q×s array, time t innovations. +! prhs[3] ghx [double] m×n array, first order reduced form. +! prhs[4] ghu [double] m×q array, first order reduced form. +! prhs[5] ghxx [double] m×n² array, second order reduced form. +! prhs[6] ghuu [double] m×q² array, second order reduced form. +! prhs[7] ghxu [double] m×nq array, second order reduced form. +! prhs[8] ghs2 [double] m×1 array, second order reduced form. +! prhs[9] ghxxx [double] m×n array, third order reduced form. +! prhs[10] ghuuu [double] m×q array, third order reduced form. +! prhs[11] ghxxu [double] m×n²q array, third order reduced form. +! prhs[12] ghxuu [double] m×nq² array, third order reduced form. +! prhs[13] ghxss [double] m×n array, third order reduced form. +! prhs[14] ghuss [double] m×q array, third order reduced form. +! prhs[15] ss [double] m×1 array, deterministic steady state +! prhs[16] numthreads [double] num of threads +! prhs[17] pruning [double] pruning option ! ! Output: -! plhs[1] y3 [double] m×s array, time t+1 particles. -! plhs[2] y2 [double] 2m×s array, time t+1 particles for the pruning latent variables up to 2nd order. The first m rows concern the pruning first-order latent variables, while the last m rows concern the pruning 2nd-order latent variables +! plhs[1] y3 [double] m×s array, time t+1 particles. +! plhs[2] ylat [double] 3m×s array, time t+1 particles for the +! pruning latent variables up to the 3rd order. Rows 1 to m contain the pruned +! first order. Rows m+1 to 2*m contain the pruned second order. Rows 2*m+1 +! to 3*m contain the pruned third order. subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') use iso_c_binding use matlab_mex @@ -347,27 +369,29 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') integer(c_int), intent(in), value :: nlhs, nrhs integer :: n, m, s, q, numthreads real(real64), pointer, contiguous :: ghx(:,:), ghu(:,:), ghxx(:,:), & - &ghuu(:,:), ghxu(:,:), ghxxx(:,:), ghuuu(:,:), ghxxu(:,:), & - &ghxuu(:,:), ghxss(:,:), ghuss(:,:), yhatlat(:,:), ylat(:,:) + &ghuu(:,:), ghxu(:,:), ghxxx(:,:), ghuuu(:,:), ghxxu(:,:), & + &ghxuu(:,:), ghxss(:,:), ghuss(:,:), yhatlat(:,:), ylat(:,:) integer :: i, j, k, xx_size, uu_size, xxx_size, uuu_size, rc character(kind=c_char, len=10) :: arg_nber type(pascal_triangle) :: p - integer, allocatable :: xx_nbeq(:), xxx_nbeq(:), & - &uu_nbeq(:), uuu_nbeq(:), xx_off(:), uu_off(:), & - &xxx_off(:), uuu_off(:) + integer, allocatable :: xxx_nbeq(:), & + &uu_nbeq(:), uuu_nbeq(:), xx_off(:), uu_off(:), & + &xxx_off(:), uuu_off(:) + integer, allocatable, target :: xx_nbeq(:) type(c_pthread_t), allocatable :: threads(:) integer, allocatable, target :: routines(:) + logical :: pruning ! 0. Checking the consistency and validity of input arguments - if (nrhs /= 15 .and. nrhs /= 17) then - call mexErrMsgTxt("Must have exactly 15 inputs or 18 inputs") + if (nrhs /= 17) then + call mexErrMsgTxt("Must have exactly 17 inputs") end if if (nlhs > 2) then call mexErrMsgTxt("Too many output arguments.") end if - do i=1, max(14, nrhs-1) + do i=1,15 if (.not. (c_associated(prhs(i)) .and. mxIsDouble(prhs(i)) .and. & (.not. mxIsComplex(prhs(i))) .and. (.not. mxIsSparse(prhs(i))))) then write (arg_nber,"(i2)") i @@ -375,20 +399,26 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') end if end do - i = max(15,nrhs) - if (.not. (c_associated(prhs(i)) .and. mxIsScalar(prhs(i)) .and. & - mxIsNumeric(prhs(i)))) then - write (arg_nber,"(i2)") i - call mexErrMsgTxt("Argument " // trim(arg_nber) // " should be a numeric scalar") + if (.not. (c_associated(prhs(16)) .and. mxIsScalar(prhs(16)) .and. & + mxIsNumeric(prhs(16)))) then + call mexErrMsgTxt("Argument 16 should be a numeric scalar") end if - numthreads = int(mxGetScalar(prhs(i))) + numthreads = int(mxGetScalar(prhs(16))) if (numthreads <= 0) then - write (arg_nber,"(i2)") i - call mexErrMsgTxt("Argument " // trim(arg_nber) // " should be a positive integer") + call mexErrMsgTxt("Argument 16 should be a positive integer") end if td3%numthreads = numthreads - n = int(mxGetM(prhs(1))) ! Number of states. + if (.not. (c_associated(prhs(17)) .and. mxIsLogicalScalar(prhs(17)))) then + call mexErrMsgTxt("Argument 17 should be a logical scalar") + end if + pruning = mxGetScalar(prhs(17)) == 1._c_double + + if (pruning) then + n = int(mxGetM(prhs(1)))/3 ! Number of states. + else + n = int(mxGetM(prhs(1))) ! Number of states. + end if s = int(mxGetN(prhs(1))) ! Number of particles. q = int(mxGetM(prhs(2))) ! Number of innovations. m = int(mxGetM(prhs(3))) ! Number of elements in the union of states and observed variables. @@ -401,13 +431,13 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') &.or. (n /= mxGetN(prhs(3))) & ! Number of columns for ghx &.or. (m /= mxGetM(prhs(4))) & ! Number of rows for ghu &.or. (q /= mxGetN(prhs(4))) & ! Number of columns for ghu - &.or. (m /= mxGetM(prhs(5))) & ! Number of rows for 2nd order constant correction + deterministic steady state - &.or. (m /= mxGetM(prhs(6))) & ! Number of rows for ghxx - &.or. (n*n /= mxGetN(prhs(6))) & ! Number of columns for ghxx - &.or. (m /= mxGetM(prhs(7))) & ! Number of rows for ghuu - &.or. (q*q /= mxGetN(prhs(7))) & ! Number of columns for ghuu - &.or. (m /= mxGetM(prhs(8))) & ! Number of rows for ghxu - &.or. (n*q /= mxGetN(prhs(8))) & ! Number of columns for ghxu + &.or. (m /= mxGetM(prhs(5))) & ! Number of rows for ghxx + &.or. (n*n /= mxGetN(prhs(5))) & ! Number of columns for ghxx + &.or. (m /= mxGetM(prhs(6))) & ! Number of rows for ghuu + &.or. (q*q /= mxGetN(prhs(6))) & ! Number of columns for ghuu + &.or. (m /= mxGetM(prhs(7))) & ! Number of rows for ghxu + &.or. (n*q /= mxGetN(prhs(7))) & ! Number of columns for ghxu + &.or. (m /= mxGetM(prhs(8))) & ! Number of rows for ghs2 &.or. (m /= mxGetM(prhs(9))) & ! Number of rows for ghxxx &.or. (n*n*n /= mxGetN(prhs(9))) & ! Number of columns for ghxxx &.or. (m /= mxGetM(prhs(10))) & ! Number of rows for ghuuu @@ -420,11 +450,8 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') &.or. (n /= mxGetN(prhs(13))) & ! Number of columns for ghxss &.or. (m /= mxGetM(prhs(14))) & ! Number of rows for ghuss &.or. (q /= mxGetN(prhs(14))) & ! Number of columns for ghuss - &.or. ((nrhs == 17) & ! With pruning optional inputs - &.and. ((2*n /= mxGetM(prhs(15))) & ! Number of rows for yhat2 - &.or. (s /= mxGetN(prhs(15))) & ! Number of columns for yhat2 - &.or. (m /= mxGetM(prhs(16)))))) then ! Number of rows for ss - ! &) then + &.or. (m /= mxGetM(prhs(15))) & ! Number of rows for ss + &) then call mexErrMsgTxt("Input dimension mismatch") end if @@ -467,27 +494,30 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') call folded_offset_loop(td3%uuu_idcs, uuu_nbeq, & &uuu_off, q, 3, p) - ! 1. Storing the relevant input variables in Fortran - td3%yhat3(1:n,1:s) => mxGetPr(prhs(1)) + ! 1. Storing the relevant input variables in Fortran + if (pruning) then + yhatlat(1:3*n,1:s) => mxGetPr(prhs(1)) + td3%yhat1 => yhatlat(1:n,1:s) + td3%yhat2 => yhatlat(n+1:2*n,1:s) + td3%yhat3 => yhatlat(2*n+1:3*n,1:s) + ! td3%xx_nbeq => xx_nbeq + else + td3%yhat3(1:n,1:s) => mxGetPr(prhs(1)) + end if td3%e(1:q,1:s) => mxGetPr(prhs(2)) ghx(1:m,1:n) => mxGetPr(prhs(3)) ghu(1:m,1:q) => mxGetPr(prhs(4)) - td3%constant => mxGetPr(prhs(5)) - ghxx(1:m,1:n*n) => mxGetPr(prhs(6)) - ghuu(1:m,1:q*q) => mxGetPr(prhs(7)) - ghxu(1:m,1:n*q) => mxGetPr(prhs(8)) + ghxx(1:m,1:n*n) => mxGetPr(prhs(5)) + ghuu(1:m,1:q*q) => mxGetPr(prhs(6)) + ghxu(1:m,1:n*q) => mxGetPr(prhs(7)) + td3%ghs2 => mxGetPr(prhs(8)) ghxxx(1:m,1:n*n*n) => mxGetPr(prhs(9)) ghuuu(1:m,1:q*q*q) => mxGetPr(prhs(10)) ghxxu(1:m,1:n*n*q) => mxGetPr(prhs(11)) ghxuu(1:m,1:n*q*q) => mxGetPr(prhs(12)) ghxss(1:m,1:n) => mxGetPr(prhs(13)) ghuss(1:m,1:q) => mxGetPr(prhs(14)) - if (nrhs == 17) then - yhatlat(1:2*n,1:s) => mxGetPr(prhs(15)) - td3%yhat1 => yhatlat(1:n,1:s) - td3%yhat2 => yhatlat(n+1:2*n,1:s) - td3%ss => mxGetPr(prhs(16)) - end if + td3%ss => mxGetPr(prhs(15)) ! Getting a transposed folded copy of the unfolded tensors ! for future loops to be more efficient @@ -543,25 +573,26 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction') plhs(1) = mxCreateDoubleMatrix(int(m, mwSize), int(s, mwSize), mxREAL) td3%y3(1:m,1:s) => mxGetPr(plhs(1)) - if (nrhs == 17) then - plhs(2) = mxCreateDoubleMatrix(int(2*m, mwSize), int(s, mwSize), mxREAL) - ylat(1:2*m,1:s) => mxGetPr(plhs(2)) - td3%y1 => ylat(1:m,1:s) - td3%y2 => ylat(m+1:2*m,1:s) + if (pruning) then + plhs(2) = mxCreateDoubleMatrix(int(3*m, mwSize), int(s, mwSize), mxREAL) + ylat(1:3*m,1:s) => mxGetPr(plhs(2)) + td3%ylat1 => ylat(1:m,1:s) + td3%ylat2 => ylat(m+1:2*m,1:s) + td3%ylat3 => ylat(2*m+1:3*m,1:s) end if allocate(threads(numthreads), routines(numthreads)) routines = [ (i, i = 1, numthreads) ] if (numthreads == 1) then - if (nrhs == 17) then + if (pruning) then call thread_eval_3_pruning(c_loc(routines(1))) else call thread_eval_3(c_loc(routines(1))) end if else ! Creating the threads - if (nrhs == 17) then + if (pruning) then do i = 1, numthreads rc = c_pthread_create(threads(i), c_null_ptr, c_funloc(thread_eval_3_pruning), c_loc(routines(i))) end do diff --git a/tests/Makefile.am b/tests/Makefile.am index 9d3704d09..2d974779e 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -35,6 +35,7 @@ MODFILES = \ analytic_derivatives/burnside_3_order_PertParamsDerivs.mod \ pruning/AnSchorfheide_pruned_state_space.mod \ pruning/AS_pruned_state_space_red_shock.mod \ + pruning/fs2000_pruning.mod \ measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod \ TeX/fs2000_corr_ME.mod \ estimation/MH_recover/fs2000_recover_tarb.mod \ diff --git a/tests/particle/first_spec_MCMC.mod b/tests/particle/first_spec_MCMC.mod index 9d6bdda20..0bae888b5 100644 --- a/tests/particle/first_spec_MCMC.mod +++ b/tests/particle/first_spec_MCMC.mod @@ -16,7 +16,7 @@ stoch_simul(order=3,periods=200, irf=0); save('my_data_MCMC.mat','ca','b'); options_.pruning=1; -options_.particle.pruning=1; +options_.particle.pruning=true; options_.particle.number_of_particles=500; estimation(datafile='my_data_MCMC.mat',order=2,mh_replic=100,filter_algorithm=sis,nonlinear_filter_initialization=2 diff --git a/tests/particle/linear_model.mod b/tests/particle/linear_model.mod index b14e2bc28..acaab56d8 100644 --- a/tests/particle/linear_model.mod +++ b/tests/particle/linear_model.mod @@ -57,7 +57,7 @@ varobs y, z; options_.particle.status = 1; options_.particle.initialization = 1; -options_.particle.pruning = 0; +options_.particle.pruning = false; options_.particle.number_of_particles = 20000; options_.particle.resampling.status = 'systematic'; options_.particle.resampling.neff_threshold = .1; diff --git a/tests/particle/local_state_space_iteration_3_test.mod b/tests/particle/local_state_space_iteration_3_test.mod index 1d1c9298f..287a3db87 100644 --- a/tests/particle/local_state_space_iteration_3_test.mod +++ b/tests/particle/local_state_space_iteration_3_test.mod @@ -40,11 +40,14 @@ dr = oo_.dr; // “rf” stands for “Reduced Form” rf_ghx = dr.ghx(dr.restrict_var_list, :); rf_ghu = dr.ghu(dr.restrict_var_list, :); +rf_ss = dr.ys(dr.order_var); +rf_ss = rf_ss(dr.restrict_var_list, :); rf_constant = dr.ys(dr.order_var)+0.5*dr.ghs2; rf_constant = rf_constant(dr.restrict_var_list, :); rf_ghxx = dr.ghxx(dr.restrict_var_list, :); rf_ghuu = dr.ghuu(dr.restrict_var_list, :); rf_ghxu = dr.ghxu(dr.restrict_var_list, :); +rf_ghs2 = dr.ghs2(dr.restrict_var_list, :); rf_ghxxx = dr.ghxxx(dr.restrict_var_list, :); rf_ghuuu = dr.ghuuu(dr.restrict_var_list, :); rf_ghxxu = dr.ghxxu(dr.restrict_var_list, :); @@ -53,7 +56,7 @@ rf_ghxss = dr.ghxss(dr.restrict_var_list, :); rf_ghuss = dr.ghuss(dr.restrict_var_list, :); % Without pruning -tStart1 = tic; for i=1:nsims, ynext1 = local_state_space_iteration_3(yhat, epsilon, rf_ghx, rf_ghu, rf_constant, rf_ghxx, rf_ghuu, rf_ghxu, rf_ghxxx, rf_ghuuu, rf_ghxxu, rf_ghxuu, rf_ghxss, rf_ghuss, options_.threads.local_state_space_iteration_3); end, tElapsed1 = toc(tStart1); +tStart1 = tic; for i=1:nsims, ynext1 = local_state_space_iteration_3(yhat, epsilon, rf_ghx, rf_ghu, rf_ghxx, rf_ghuu, rf_ghxu, rf_ghs2, rf_ghxxx, rf_ghuuu, rf_ghxxu, rf_ghxuu, rf_ghxss, rf_ghuss, rf_ss, options_.threads.local_state_space_iteration_3, false); end, tElapsed1 = toc(tStart1); tStart2 = tic; [udr] = folded_to_unfolded_dr(dr, M_, options_); for i=1:nsims, ynext2 = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_, udr); end, tElapsed2 = toc(tStart2); if max(max(abs(ynext1-ynext2))) > 1e-10 @@ -71,19 +74,20 @@ else end % With pruning -rf_ss = dr.ys(dr.order_var); -rf_ss = rf_ss(dr.restrict_var_list, :); -yhat_ = chol(oo_.var(state_idx,state_idx))*randn(nstates, nparticles); -yhat__ = zeros(2*nstates,nparticles); -yhat__(1:nstates,:) = yhat_; -yhat__(nstates+1:2*nstates,:) = chol(oo_.var(state_idx,state_idx))*randn(nstates, nparticles); +yhat1_KKSS = chol(oo_.var(state_idx,state_idx))*randn(nstates, nparticles); +yhat2_KKSS = chol(oo_.var(state_idx,state_idx))*randn(nstates, nparticles); +yhat1_AVR = yhat1_KKSS; +yhat2_AVR = yhat2_KKSS-yhat1_AVR; +yhat3_AVR = chol(oo_.var(state_idx,state_idx))*randn(nstates, nparticles); +yhat_AVR = [yhat1_AVR;yhat2_AVR;yhat3_AVR]; nstatesandobs = size(rf_ghx,1); -[ynext1,ynext1_lat] = local_state_space_iteration_3(yhat, epsilon, rf_ghx, rf_ghu, rf_constant, rf_ghxx, rf_ghuu, rf_ghxu, rf_ghxxx, rf_ghuuu, rf_ghxxu, rf_ghxuu, rf_ghxss, rf_ghuss, yhat__, rf_ss, options_.threads.local_state_space_iteration_3); -[ynext2,ynext2_lat] = local_state_space_iteration_2(yhat__(nstates+1:2*nstates,:), epsilon, rf_ghx, rf_ghu, rf_constant, rf_ghxx, rf_ghuu, rf_ghxu, yhat_, rf_ss, options_.threads.local_state_space_iteration_2); -if max(max(abs(ynext1_lat(nstatesandobs+1:2*nstatesandobs,:)-ynext2))) > 1e-14 +[ynext3,ynext3_lat] = local_state_space_iteration_3(yhat_AVR, epsilon, rf_ghx, rf_ghu, rf_ghxx, rf_ghuu, rf_ghxu, rf_ghs2, rf_ghxxx, rf_ghuuu, rf_ghxxu, rf_ghxuu, rf_ghxss, rf_ghuss, rf_ss, options_.threads.local_state_space_iteration_3, true); +[ynext2,ynext2_lat] = local_state_space_iteration_2(yhat2_KKSS, epsilon, rf_ghx, rf_ghu, rf_constant, rf_ghxx, rf_ghuu, rf_ghxu, yhat1_KKSS, rf_ss, options_.threads.local_state_space_iteration_2); +if max(max(abs(ynext3_lat(1:nstatesandobs,:)-ynext2_lat))) > 1e-14 + error('With pruning: inconsistency between local_state_space_iteration_2 and local_state_space_iteration_3 on the 1st-order pruned variable') +end +ynext3_2nd = ynext3_lat(1:nstatesandobs,:)+ynext3_lat(nstatesandobs+1:2*nstatesandobs,:)-rf_ss; +if max(abs(ynext3_2nd(:)-ynext2(:))) > 1e-14 error('With pruning: inconsistency between local_state_space_iteration_2 and local_state_space_iteration_3 on the 2nd-order pruned variable') end -if max(max(abs(ynext1_lat(1:nstatesandobs,:)-ynext2_lat))) > 1e-14 - error('With pruning: inconsistency between local_state_space_iteration_2 and local_state_space_iteration_3 on the 1st-order pruned variable') -end \ No newline at end of file diff --git a/tests/particle/noisy_ar1.mod b/tests/particle/noisy_ar1.mod index aaebcfea2..798f1b881 100644 --- a/tests/particle/noisy_ar1.mod +++ b/tests/particle/noisy_ar1.mod @@ -39,7 +39,7 @@ particle = 1; if particle; options_.particle.status = 1; options_.particle.initialization = 1; -options_.particle.pruning = 0; +options_.particle.pruning = false; options_.particle.number_of_particles = 20000; options_.particle.resampling.status = 'systematic'; options_.particle.resampling.neff_threshold = .1; diff --git a/tests/pruning/fs2000_pruning.mod b/tests/pruning/fs2000_pruning.mod new file mode 100644 index 000000000..e39a42606 --- /dev/null +++ b/tests/pruning/fs2000_pruning.mod @@ -0,0 +1,198 @@ +/* +Compare simulation results with pruning at second order + */ + +var m P c e W R k d n l gy_obs gp_obs y dA ; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +steady_state_model; + dA = exp(gam); + gst = 1/dA; + m = mst; + khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1)); + xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1); + nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp ); + n = xist/(nust+xist); + P = xist + nust; + k = khst*n; + + l = psi*mst*n/( (1-psi)*(1-n) ); + c = mst/P; + d = l - mst + 1; + y = k^alp*n^(1-alp)*gst^alp; + R = mst/bet; + W = l/n; + ist = y-c; + q = 1 - d; + + e = 1; + + gp_obs = m/dA; + gy_obs = dA; +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +T = 10; +ex=randn(T,M_.exo_nbr); + +% 2nd order +stoch_simul(order=2, nograph, irf=0); + +% simult_.m: implements KKSS +options_.pruning = true; +Y2_simult = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); + +% local_state_space_iteration_2 mex: implements KKSS +constant = oo_.dr.ys(oo_.dr.order_var)+0.5*oo_.dr.ghs2; +ss = oo_.dr.ys(oo_.dr.order_var); +Y1_local = zeros(M_.endo_nbr,T+1); +Y1_local(:,1) = oo_.dr.ys; +Y2_local = zeros(M_.endo_nbr,T+1); +Y2_local(:,1) = oo_.dr.ys; +state_var = oo_.dr.order_var(M_.nstatic+1:M_.nstatic+M_.nspred); +for t=2:T+1 + u = ex(t-1,:)'; + yhat1 = Y1_local(state_var,t-1)-oo_.dr.ys(state_var); + yhat2 = Y2_local(state_var,t-1)-oo_.dr.ys(state_var); + [Y2, Y1] = local_state_space_iteration_2(yhat2, u, oo_.dr.ghx, oo_.dr.ghu, constant, oo_.dr.ghxx, oo_.dr.ghuu, oo_.dr.ghxu, yhat1, ss, options_.threads.local_state_space_iteration_2); + Y1_local(oo_.dr.order_var,t) = Y1; + Y2_local(oo_.dr.order_var,t) = Y2; +end + +if (max(abs(Y2_local(:) - Y2_simult(:)))>1e-12) + error("2nd-order output of simult_ and local_state_space_iteration_2 output are inconsistent.") +end + +% pruned_state_space_system.m: implements Andreasen et al. +pss = pruned_state_space_system(M_, options_, oo_.dr, [], 0, false, false); +Y2_an = zeros(M_.endo_nbr,T+1); +Y2_an(:,1) = oo_.dr.ys; +% z = [xf;xs;kron(xf,xf)] +% inov = [u;kron(u,u)-E_uu(:);kron(xf,u)] +z = zeros(2*M_.nspred+M_.nspred^2, 1); +z(1:2*M_.nspred, 1) = repmat(Y2_an(state_var,1)-oo_.dr.ys(state_var), 2, 1); +z(2*M_.nspred+1:end, 1) = kron(z(1:M_.nspred, 1), z(M_.nspred+1:2*M_.nspred, 1)); +inov = zeros(M_.exo_nbr+M_.exo_nbr^2+M_.nspred*M_.exo_nbr, 1); +for t=2:T+1 + u = ex(t-1,:)'; + inov(1:M_.exo_nbr, 1) = u; + inov(M_.exo_nbr+1:M_.exo_nbr+M_.exo_nbr^2, 1) = kron(u, u) - M_.Sigma_e(:); + inov(M_.exo_nbr+M_.exo_nbr^2+1:end, 1) = kron(z(1:M_.nspred, 1), u); + Y = oo_.dr.ys(oo_.dr.order_var) + pss.d + pss.C*z + pss.D*inov; + x1 = z(1:M_.nspred,1); + x2 = z(M_.nspred+1:2*M_.nspred,1); + z = pss.c + pss.A*z + pss.B*inov; + Y2_an(oo_.dr.order_var,t) = Y; +end + +if (max(abs(Y2_an(:) - Y2_local(:)))>1e-12) + error("2nd-order output of pruned_state_space_system and local_state_space_iteration_2 output are inconsistent.") +end + +% 3rd order +stoch_simul(order=3, nograph, irf=0); +% simult_.m +Y3_simult = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); +% pruned_state_space_system.m +pss = pruned_state_space_system(M_, options_, oo_.dr, [], 0, false, false); +Y3_an = zeros(M_.endo_nbr,T+1); +Y3_an(:,1) = oo_.dr.ys; +% z = [xf; xs; kron(xf,xf); xrd; kron(xf,xs); kron(kron(xf,xf),xf)] +% inov = [u; kron(u,u)-E_uu(:); kron(xf,u); kron(xs,u); kron(kron(xf,xf),u); kron(kron(xf,u),u); kron(kron(u,u),u))] +x_nbr = M_.nspred; +u_nbr = M_.exo_nbr; +z_nbr = x_nbr + x_nbr + x_nbr^2 + x_nbr + x_nbr^2 + x_nbr^3; +inov_nbr = u_nbr + u_nbr^2 + x_nbr*u_nbr + x_nbr*u_nbr + x_nbr^2*u_nbr + x_nbr*u_nbr^2 + u_nbr^3; +z = zeros(z_nbr, 1); +inov = zeros(inov_nbr, 1); +for t=2:T+1 + u = ex(t-1,:)'; + xf = z(1:M_.nspred,1); + xs = z(M_.nspred+1:2*M_.nspred,1); + inov(1:u_nbr, 1) = u; + ub = u_nbr; + inov(ub+1:ub+u_nbr^2, 1) = kron(u, u) - M_.Sigma_e(:); + ub = ub+u_nbr^2; + inov(ub+1:ub+x_nbr*u_nbr, 1) = kron(xf, u); + ub = ub+x_nbr*u_nbr; + inov(ub+1:ub+x_nbr*u_nbr, 1) = kron(xs, u); + ub = ub+x_nbr*u_nbr; + inov(ub+1:ub+x_nbr^2*u_nbr) = kron(kron(xf,xf),u); + ub = ub+x_nbr^2*u_nbr; + inov(ub+1:ub+x_nbr*u_nbr^2) = kron(kron(xf,u),u); + ub = ub+x_nbr*u_nbr^2; + inov(ub+1:end) = kron(kron(u,u),u); + Y = oo_.dr.ys(oo_.dr.order_var) + pss.d + pss.C*z + pss.D*inov; + z = pss.c + pss.A*z + pss.B*inov; + Y3_an(oo_.dr.order_var,t) = Y; +end + +if (max(abs(Y3_an(:) - Y3_simult(:)))>1e-12) + error("3rd-order outputs of pruned_state_space_system and simult_ are inconsistent.") +end + +% local_state_space_iteration_3 mex +Y3_local = zeros(M_.endo_nbr,T+1); +Y3_local(:,1) = oo_.dr.ys; +Ylat_local = zeros(3*M_.endo_nbr,T+1); +Ylat_local(:,1) = repmat(oo_.dr.ys,3,1); +for t=2:T+1 + u = ex(t-1,:)'; + yhat1 = Ylat_local(state_var,t-1)-oo_.dr.ys(state_var); + yhat2 = Ylat_local(M_.endo_nbr+state_var,t-1)-oo_.dr.ys(state_var); + yhat3 = Ylat_local(2*M_.endo_nbr+state_var,t-1)-oo_.dr.ys(state_var); + ylat = [yhat1;yhat2;yhat3]; + [Y3,Y] = local_state_space_iteration_3(ylat, u, oo_.dr.ghx, oo_.dr.ghu, oo_.dr.ghxx, oo_.dr.ghuu, oo_.dr.ghxu, oo_.dr.ghs2, oo_.dr.ghxxx, oo_.dr.ghuuu, oo_.dr.ghxxu, oo_.dr.ghxuu, oo_.dr.ghxss, oo_.dr.ghuss, ss, options_.threads.local_state_space_iteration_3, true); + Ylat_local(oo_.dr.order_var,t) = Y(1:M_.endo_nbr); + Ylat_local(M_.endo_nbr+oo_.dr.order_var,t) = Y(M_.endo_nbr+1:2*M_.endo_nbr); + Ylat_local(2*M_.endo_nbr+oo_.dr.order_var,t) = Y(2*M_.endo_nbr+1:end); + Y3_local(oo_.dr.order_var,t) = Y3; +end + +Y3_local_1 = Ylat_local(1:M_.endo_nbr,:); +Y3_local_2 = Y3_local_1 + Ylat_local(M_.endo_nbr+1:2*M_.endo_nbr,:) - oo_.dr.ys; + +if (max(abs(Y3_local_1(:) - Y1_local(:)))>1e-12) + error("1st-order outputs of local_state_space_iteration_3 and local_state_space_iteration_2 are inconsistent.") +end + +if (max(abs(Y3_local_2(:) - Y2_local(:)))>4e-12) + error("2nd-order outputs of local_state_space_iteration_3 and local_state_space_iteration_2 are inconsistent.") +end + +if (max(abs(Y3_local(:) - Y3_simult(:)))>1e-12) + error("3rd-order output of simult_ and local_state_space_iteration_3 output are inconsistent.") +end