Improve debug mode.

meson.build

@@ -353,6 +353,8 @@ shared_module('qmc_sequence', qmc_sequence_src, kwargs : qmc_sequence_mex_kwargs
 
 shared_module('kalman_steady_state', 'mex/sources/kalman_steady_state/kalman_steady_state.cc', kwargs : mex_kwargs, dependencies : slicot_dep)
 
+shared_module('kalman_filter', [ 'mex/sources/kalman_filter/kalman_filter.f08' ] + mex_blas_fortran_iface,
+              kwargs : mex_kwargs, dependencies : [ blas_dep, lapack_dep ])
 
 ## k-order simulation stuff
 
@@ -1752,6 +1754,7 @@ mod_and_m_tests = [
   { 'test' : [ 'cyclereduction.m' ] },
   { 'test' : [ 'logarithmicreduction.m' ] },
   { 'test' : [ 'riccatiupdate.m' ] },
+  { 'test' : [ 'kalman/likelihood/test_kalman_mex.m' ] },
   { 'test' : [ 'contribs.m' ],
     'extra' : [ 'sandbox.mod',
                 'simulateddata.m' ] },

mex/sources/blas_lapack.F08

@@ -145,6 +145,38 @@ module lapack
        integer(blint), intent(out) :: sdim, info
      end subroutine dgges
 
+     subroutine dgecon(norm, n, a, lda, anorm, rcond, work, iwork, info)
+       import :: blint, real64
+       implicit none
+       character, intent(in) :: norm
+       integer(blint), intent(in) :: n, lda
+       real(real64), dimension(*), intent(in) :: a
+       real(real64), intent(in) :: anorm
+       real(real64), intent(out) :: rcond
+       real(real64), dimension(*), intent(out) :: work
+       integer(blint), dimension(*), intent(out) :: iwork
+       integer(blint), intent(out) :: info
+     end subroutine dgecon
+
+     subroutine dgetrf(m, n, a, lda, ipiv, info)
+       import :: blint, real64
+       implicit none
+       integer(blint), intent(in) :: m, n, lda
+       real(real64), dimension(*), intent(inout) :: a
+       integer(blint), dimension(*), intent(out) :: ipiv
+       integer(blint), intent(out) :: info
+     end subroutine dgetrf
+
+     subroutine dgetri(n, a, lda, ipiv, work, lwork, info)
+       import :: blint, real64
+       implicit none
+       integer(blint), intent(in) :: n, lda, lwork
+       real(real64), dimension(*), intent(inout) :: a
+       integer(blint), dimension(*), intent(in):: ipiv
+       real(real64), dimension(*), intent(out) :: work
+       integer(blint), intent(out) :: info
+     end subroutine dgetri
+
      subroutine dpotrf(uplo, n, a, lda, info)
        import :: blint, real64
        implicit none
@@ -163,6 +195,14 @@ module lapack
        real(real64), dimension(*), intent(out) :: work
        real(real64) :: dlange
      end function dlange
+
+     function ilaenv(ispec, name, opts, n1, n2, n3, n4)
+       import :: blint  
+       implicit none
+       integer(blint) :: ispec, n1, n2, n3, n4
+       character(len=*), intent(in) :: name, opts
+       integer(blint) :: ilaenv
+     end function ilaenv
   end interface
 
 contains

mex/sources/kalman_filter/kalman_filter.f08

@@ -0,0 +1,435 @@
+! Copyright © 2023 Dynare Team
+!
+! This file is part of Dynare.
+!
+! Dynare is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! Dynare is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+!
+! input:
+! prhs[1]  Y                  [double] p x T data matrix
+! prhs[2]  start              [integer] first period index
+! prhs[3]  last               [integer] last period index
+! prhs[4]  a                  [double] m x 1 initial mean of the state vector
+! prhs[5]  P                  [double] m x m initial variance-covariance matrix
+!                                            of the state vector
+! prhs[6]  kalman_tol         [double] tolerance parameter (rcond, inversibility
+!                                     of the covariance matrix of the prediction
+!                                     errors)
+! prhs[7]  riccati_tol        [double] tolerance parameter (iteration over the
+!                                      Riccati equation). 
+! prhs[8]  presample          [integer] presampling if strictly positive (number
+!                                      of initial iterations to be discarded
+!                                      when evaluating the likelihood).
+! prhs[9]  T                  [double] m x m transition matrix of the state
+!                                      equation 
+! prhs[10] Q                  [double] r x r variance-covariance matrix of the
+!                                      structural innovations (noise in the 
+!                                      state equation).
+! prhs[11] R                  [double] m x r second matrix of the state equation
+!                                      relating the structural innovations to
+!                                      the state variables.
+! prhs[12] H                  [double] p x p variance-covariance matrix of the
+!                                      measurement errors (if no measurement
+!                                      errors set H as a zero scalar).
+! prhs[13] Z                  [double] p x m matrix relating the states to the
+!                                      observed variables or vector of indices
+!                                      (depending on the value of Zflag).
+! prhs[14] Zflag             [integer] equal to 0 if Z is a vector of indices
+!                                      targeting the observed variables in the
+!                                      state vector, equal to 1 if Z is a
+!                                      p x m matrix.
+! prhs[15] diffuse_periods   [integer] number of diffuse filter
+!                                      periods in the initialization step.   
+! output:
+! plhs[1]  LIK                [double] value of (minus) the likelihood.
+! plhs[2]  LIKK               [double] vector containing the density of each
+!                                      observation.
+! plhs[3]  a                  [double] m x 1 mean of the state vector at the end
+!                                      of the (sub)sample.
+! plhs[4]  P                  [double] m x m variance-covariance matrix of the
+!                                      state vector at the end of the
+!                                      (sub)sample.
+! The current Kalman filter is called by dsge_likelihood.m when
+! - kalman_algo = 1 or 3 (multivariate Kalman filter)
+! - no missing data
+! - no use of occbin
+! - diffuse_periods = False
+! - fast_kalman_filter = True
+
+subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name="mexFunction")
+   use matlab_mex
+   use blas
+   use lapack
+   use iso_c_binding
+   implicit none
+
+   type(c_ptr), dimension(*), intent(in), target :: prhs
+   type(c_ptr), dimension(*), intent(out) :: plhs
+   integer(c_int), intent(in), value :: nlhs, nrhs
+
+   type(c_ptr) :: Y_mx, a_mx, P_mx, kalman_tol_mx, &
+                 &riccati_tol_mx, presample_mx, T_mx, Q_mx, R_mx, H_mx, &
+                 &Z_mx, Zflag_mx, diffuse_periods_mx
+   type(c_ptr), dimension(2) :: call_lhs
+   type(c_ptr), dimension(11) :: call_rhs
+   integer :: presample, diffuse_periods, t, nper, smpl, s
+   integer, allocatable, dimension(:) :: indZ
+   integer(c_int) :: retval
+   integer(blint) :: m, r, p, info, lwork, nb, i, j
+   real(real64) :: kalman_tol, riccati_tol, rcond, log_dF, lik, pi 
+   real(real64), dimension(:,:), contiguous, pointer :: PP, TT, Q, RR, H, Z, Y, K, K_m, iF_m
+   real(real64), dimension(:), contiguous, pointer :: a, likk, a_m
+   logical :: Zflag, steady_flag, Hflag
+   real(real64), allocatable :: v(:), F(:,:), lu(:,:), QQ(:,:),          &
+                               &work_rcond(:), work_inv(:), P_next(:,:), &
+                               &tmp_m_m(:,:), tmp_m_p(:,:), tmp_m_r(:,:),&
+                               &tmp_m_m_prime(:,:), tmp_v(:), a_next(:), &
+                               &old_K(:,:), P_iter(:,:), a_iter(:)
+   integer(blint), allocatable :: iwork(:), ipiv(:) 
+
+   ! Check the number of output arguments
+   if (nlhs < 2) then
+      call mexErrMsgTxt("Must have at least 2 outputs")
+   end if
+
+   ! Check the consistency and validity of input arguments
+   if (nrhs < 9) then
+      call mexErrMsgTxt("Must have at least 9 inputs")
+   end if
+
+   Y_mx = prhs(1)
+   a_mx = prhs(2)
+   P_mx = prhs(3)
+   kalman_tol_mx = prhs(4)
+   riccati_tol_mx = prhs(5)
+   T_mx = prhs(6)
+   Q_mx = prhs(7)
+   R_mx = prhs(8)
+   Z_mx = prhs(9)
+
+   if (.not. mxIsDouble(Y_mx) .or. mxIsComplex(Y_mx) .or. mxIsSparse(Y_mx)) then
+      call mexErrMsgTxt("1st argument (Y) should be a real dense matrix")
+   end if
+   if (.not. (mxIsDouble(a_mx) .and. ((mxGetM(a_mx) == 1) .or. &
+      &(mxGetN(a_mx) == 1))) .or. mxIsComplex(a_mx) .or. mxIsSparse(a_mx)) then
+      call mexErrMsgTxt("2nd argument (a) should be a real dense vector")
+   end if
+   if (.not. mxIsDouble(P_mx) .or. mxIsComplex(P_mx) .or. mxIsSparse(P_mx)) then
+      call mexErrMsgTxt("3rd argument (P) should be a real dense matrix")
+   end if
+   if (.not. (mxIsScalar(kalman_tol_mx) .and. mxIsNumeric(kalman_tol_mx))) then
+      call mexErrMsgTxt("4th argument (kalman_tol) should be a numeric scalar")
+   end if
+   if (.not. (mxIsScalar(riccati_tol_mx) .and. mxIsNumeric(riccati_tol_mx))) then
+      call mexErrMsgTxt("5th argument (riccati_tol) should be a numeric scalar")
+   end if
+   if (.not. mxIsDouble(T_mx) .or. mxIsComplex(T_mx) .or. mxIsSparse(T_mx)) then
+      call mexErrMsgTxt("6th argument (T) should be a real dense matrix")
+   end if
+   if (.not. mxIsDouble(Q_mx) .or. mxIsComplex(Q_mx) .or. mxIsSparse(Q_mx)) then
+      call mexErrMsgTxt("7th argument (Q) should be a real dense matrix")
+   end if
+   if (.not. mxIsDouble(R_mx) .or. mxIsComplex(R_mx) .or. mxIsSparse(R_mx)) then
+      call mexErrMsgTxt("8th argument (R) should be a real dense matrix")
+   end if
+
+   ! Import variables
+   a => mxGetPr(a_mx)
+   m = size(a, 1, blint)
+   p = int(mxGetM(Y_mx), blint)
+   nper = int(mxGetN(Y_mx))
+   r = int(mxGetM(Q_mx), blint)
+   Y(1:p,1:nper) => mxGetPr(Y_mx)
+   PP(1:m,1:m) => mxGetPr(P_mx)
+   TT(1:m,1:m) => mxGetPr(T_mx)
+   Q(1:r,1:r) => mxGetPr(Q_mx)
+   RR(1:m,1:r) => mxGetPr(R_mx)
+   kalman_tol = mxGetScalar(kalman_tol_mx)
+   riccati_tol = mxGetScalar(riccati_tol_mx)
+
+   ! Check the order consistency of input matrices
+   if ((mxGetM(P_mx) /= m) .or. (mxGetN(P_mx) /= m) .or.& ! P
+      &(mxGetM(T_mx) /= m) .or. (mxGetN(T_mx) /= m) .or.& ! T
+      &(mxGetN(Q_mx) /= r) .or.                         & ! Q
+      &(mxGetM(R_mx) /= m) .or. (mxGetN(R_mx) /= r)) then ! R
+      call mexErrMsgTxt("Input dimension mismatch in (a, Y, P, T, Q, R)")
+   end if
+
+   ! Optional inputs
+   ! Zflag
+   if (nrhs > 9) then
+      Zflag_mx = prhs(10)
+      if (.not. (mxIsScalar(Zflag_mx) .and. mxIsNumeric(Zflag_mx))) then
+         call mexErrMsgTxt("10th argument (Zflag) should be a numeric scalar")
+      end if
+      Zflag = (mxGetScalar(Zflag_mx) == 1._c_int)
+   else
+      Zflag = .false.
+   end if
+   if (Zflag) then
+      if (.not. mxIsDouble(Z_mx) .or. mxIsComplex(Z_mx) .or. mxIsSparse(Z_mx)) then
+         call mexErrMsgTxt("9th argument (Z) should be a real dense matrix")
+      end if
+      if ((mxGetM(Z_mx) /= p) .or. (mxGetN(Z_mx) /= m)) then
+         call mexErrMsgTxt("Input dimension mismatch in Z")
+      end if
+      Z(1:p,1:m) => mxGetPr(Z_mx)
+   else
+      if (.not. (mxIsDouble(Z_mx) .and. ((mxGetM(Z_mx) == 1) .or. &  
+         &(mxGetN(Z_mx) == 1))) .or. mxIsComplex(Z_mx) .or.       &
+         &mxIsSparse(Z_mx)) then
+         call mexErrMsgTxt("9th argument (Z) should be a real dense vector")
+      end if
+      if ((mxGetM(Z_mx) /= p) .and. (mxGetN(Z_mx) /= p)) then
+         call mexErrMsgTxt("Input dimension mismatch in Z")
+      end if
+      Z(1:p,1:1) => mxGetPr(Z_mx)
+      indZ = int(Z(1:p,1))
+   end if
+
+   ! H
+   if (nrhs > 10) then
+      H_mx = prhs(11)
+      if ((mxIsScalar(H_mx) .and. mxIsNumeric(H_mx))) then
+         Hflag = .false.
+      elseif (mxIsDouble(H_mx) .and. .not. (mxIsComplex(H_mx) .or. mxIsSparse(H_mx))) then
+         Hflag = .true.
+         if ((mxGetM(H_mx) /= p) .or. (mxGetN(H_mx) /= p)) then
+            call mexErrMsgTxt("Input dimension mismatch in H")
+         end if
+         H(1:p,1:p) => mxGetPr(H_mx)
+      else
+         call mexErrMsgTxt("11th argument (H) should be a real dense matrix or a zero scalar if no measurement error is set")
+      end if
+   else
+      Hflag = .false.
+   end if
+
+   ! diffuse_periods
+   if (nrhs > 11) then
+      diffuse_periods_mx = prhs(12)
+      if (.not. (mxIsScalar(diffuse_periods_mx) .and. mxIsNumeric(diffuse_periods_mx))) then
+         call mexErrMsgTxt("12th argument (diffuse_periods) should be a numeric scalar")
+      end if
+      diffuse_periods = int(mxGetScalar(diffuse_periods_mx))
+   else
+      diffuse_periods = 0
+   end if
+
+   ! presample
+   if (nrhs > 12) then
+      presample_mx = prhs(13)
+      if (.not. (mxIsScalar(presample_mx) .and. mxIsNumeric(presample_mx))) then
+         call mexErrMsgTxt("13th argument (presample) should be a numeric scalar")
+      end if
+      presample = int(mxGetScalar(presample_mx))
+   else
+      presample = 0
+   end if
+
+   smpl = nper-diffuse_periods
+
+   ! Density of each observation
+   plhs(2) = mxCreateDoubleMatrix(int(smpl, mwSize), 1_mwSize, mxREAL)
+   likk(1:smpl) => mxGetPr(plhs(2))
+
+   ! Optimal number of blocks for DGETRI
+   nb = ilaenv(1_blint, "DGETRI", " ", p, -1_blint, -1_blint, -1_blint)
+   lwork = p*nb
+
+   allocate(v(p), F(p,p), ipiv(p), work_rcond(4*p), work_inv(lwork), iwork(p), &
+           &tmp_m_m(m,m), tmp_m_r(m,r), tmp_m_p(m,p), tmp_m_m_prime(m,m), &
+           &tmp_v(p), old_K(m,p), K(m,p), a_next(m), QQ(m,m))
+
+   ! Compute RQR'
+   ! (i) tmp <- RQ
+   call matmul_add("N", "N", 1._real64, RR, Q, 0._real64, tmp_m_r)
+   ! (ii) QQ <- tmp*R'
+   call matmul_add("N", "T", 1._real64, tmp_m_r, RR, 0._real64, QQ)
+
+   t = diffuse_periods+1
+   steady_flag = .false.
+   log_dF = huge(log_dF)
+   old_K = huge(0._real64)
+   a_iter = a
+   P_iter = PP
+   do 
+      if ((t > nper) .or. (steady_flag)) exit
+      s = t-diffuse_periods
+      ! v <- Y(:,t) - Z*a
+      if (Zflag) then
+         v = Y(:,t)
+         call dgemv("N", p, m, -1._real64, Z, p, a_iter, 1_blint, 1._real64, v, 1_blint)
+         ! F <- Z*P*Z' + H
+         ! (i) tmp <- P*Z'
+         call matmul_add("N", "T", 1._real64, P_iter, Z, 0._real64, tmp_m_p)
+         ! (ii) F <- Z*tmp + H
+         if (Hflag) then
+            F = H
+            call matmul_add("N", "N", 1._real64, Z, tmp_m_p, 1._real64, F)
+         else
+            call matmul_add("N", "N", 1._real64, Z, tmp_m_p, 0._real64, F)
+         end if
+      else
+         v = Y(:,t) - a_iter(indZ)
+         if (Hflag) then
+            F = P_iter(indZ,indZ) + H
+         else
+            F = P_iter(indZ,indZ)
+         end if
+      end if
+
+      ! Compute the reciprocal condition number of F using its LU 
+      ! decomposition
+      ! (i) LU decomposition of F
+      lu = F
+      call dgetrf(p, p, lu, p, ipiv, info)
+      if (info < 0) then
+         call mexErrMsgTxt("Ill-conditioned F!")
+      end if
+      ! (ii) Reciprocal condition number of F
+      call dgecon("1", p, lu, p, norm(F, "1"), rcond, work_rcond, iwork, info)
+      if (rcond < kalman_tol) then
+         call mexErrMsgTxt("Ill-conditioned F!")
+      end if
+
+      ! Compute log(det(F))
+      log_dF = 0._real64
+      do i=1,p
+         log_dF = log_dF+log(abs(lu(i,i)))
+      end do
+      
+      ! Compute the inverse of F using its LU decomposition
+      call dgetri(p, lu, p, ipiv, work_inv, lwork, info)
+
+      ! Density of each observation
+      call dgemv("N", p, p, 1._real64, lu, p, v, 1_blint, 0._real64, tmp_v, 1_blint)
+      likk(s) = log_dF
+      do i=1,p
+         likk(s) = likk(s)+v(i)*tmp_v(i)
+      end do
+
+      ! Compute K
+      if (Zflag) then
+         ! Compute K = PZ'F^{-1}
+         ! (i) tmp <- PZ'
+         call matmul_add("N", "T", 1._real64, P_iter, Z, 0._real64, tmp_m_p)
+         ! (ii) K <- tmp*F^{-1}
+         call matmul_add("N", "N", 1._real64, tmp_m_p, lu, 0._real64, K)
+      else
+         ! Compute K = P(:,Z)F^{-1} 
+         call matmul_add("N", "N", 1._real64, P_iter(:,indZ), lu, 0._real64, K)
+      end if
+
+      P_next = QQ
+
+      ! Compute P_{t+1}
+      if (Zflag) then
+         ! Compute P_{t+1} <- T(I-KZ)PT' + RQR'
+         ! (i) Compute tmp' <- I-KZ
+         do j=1,m
+            do i=1,m
+               if (j == i) then
+                  tmp_m_m_prime(i,j) = 1._real64
+               else
+                  tmp_m_m_prime(i,j) = 0._real64
+               end if
+            end do
+         end do
+         call matmul_add("N", "N", -1._real64, K, Z, 1._real64, tmp_m_m_prime)
+         ! (ii) tmp_m_m <- tmp'*P
+         call matmul_add("N", "N", 1._real64, tmp_m_m_prime, P_iter, 0._real64, tmp_m_m)
+         ! (iii) tmp_m_m' <- T*tmp_m_m 
+         call matmul_add("N", "N", 1._real64, TT, tmp_m_m, 0._real64, tmp_m_m_prime)
+         ! (iv) P_next <- tmp_m_m'*T' + P_next
+         call matmul_add("N", "T", 1._real64, tmp_m_m_prime, TT, 1._real64, P_next)
+      else
+         ! Compute P_{t+1} <- T(P-K*P(Z,:))T' + RQR'
+         ! (i) tmp_m_m <- P-K*P(Z,:)
+         tmp_m_m = P_iter
+         call matmul_add("N", "N", -1._real64, K, P_iter(indZ,:), 1._real64, tmp_m_m)
+         ! (ii) tmp_m_m' <- T*tmp_m_m
+         call matmul_add("N", "N", 1._real64, TT, tmp_m_m, 0._real64, tmp_m_m_prime)
+         ! (iii) P_next <- tmp_m_m'*T' + P_next
+         call matmul_add("N", "T", 1._real64, tmp_m_m_prime, TT, 1._real64, P_next)
+      end if
+
+      ! Compute a_{t+1} = T*(a+Kv)
+      ! (i) a <- a+Kv
+      call dgemv("N", m, p, 1._real64, K, m, v, 1_blint, 1._real64, a_iter, 1_blint)
+      ! (ii) a_next <- T*tmp_v
+      call dgemv("N", m, m, 1._real64, TT, m, a_iter, 1_blint, 0._real64, a_next, 1_blint)
+      
+      ! Check the wedge between gain matrices
+      steady_flag = (norm(K-old_K, "M") <= riccati_tol)
+
+      ! Set up next iteration
+      old_K = K
+      a_iter = a_next
+      P_iter = P_next
+      t = t+1
+   end do
+
+   ! Add observation's densities constants and divide by two.
+   pi=4._real64*DATAN(1._real64)
+   likk(1:s) = 0.5_real64*(likk(1:s) + p*log(2*pi))
+
+   ! Call kalman_filter_ss if necessary 
+   if (t <= nper) then
+      ! call mexPrintf("Calling kalman_filter_ss!")
+      call_rhs(1) = Y_mx
+      call_rhs(2) = mxCreateDoubleScalar(real(t, c_double))
+      call_rhs(3) = mxCreateDoubleScalar(real(nper, c_double))
+      call_rhs(4) = mxCreateDoubleMatrix(int(m, mwSize), 1_mwSize, mxREAL)
+      a_m => mxGetPr(call_rhs(4))
+      a_m = a_iter
+      call_rhs(5) = T_mx
+      call_rhs(6) = mxCreateDoubleMatrix(int(m, mwSize), int(p, mwSize), mxREAL)
+      K_m(1:m,1:p) => mxGetPr(call_rhs(6)) 
+      K_m = K
+      call_rhs(7) = mxCreateDoubleMatrix(int(p, mwSize), int(p, mwSize), mxREAL)
+      iF_m(1:p,1:p) => mxGetPr(call_rhs(7)) 
+      iF_m = lu
+      call_rhs(8) = mxCreateDoubleScalar(log_dF)
+      call_rhs(9) = Z_mx
+      call_rhs(10) = mxCreateDoubleScalar(real(p, c_double))
+      if (nrhs > 9) then
+         call_rhs(11) = Zflag_mx
+      else
+         call_rhs(11) = mxCreateDoubleScalar(0._c_double)
+      end if
+      retval = mexCallMATLAB(2_c_int, call_lhs, 11_c_int, call_rhs, "kalman_filter_ss")
+      if (retval /= 0_c_int) then
+         call mexErrMsgTxt("Error calling kalman_filter_ss!")
+      end if
+      likk(s+1:smpl) = mxGetPr(call_lhs(2))
+   end if
+
+   ! Compute minus the log-likelihood.
+   if (presample > diffuse_periods) then
+      lik = sum(likk(1+presample-diffuse_periods:smpl))
+   else
+      lik = sum(likk)
+   end if
+   plhs(1) = mxCreateDoubleScalar(real(lik, c_double))
+      
+   ! N.B.: An alternative using Cholesky and the symmetry of matrices goes
+   ! (1) P = L*L' (dpotrf)
+   ! (2) F <- Z*L*(Z*L)' + H (dsyrk)
+   ! F is necessarily symmetric positive-definite in this case
+   ! (3) put F in packed format
+   ! (4) compute its Cholesky decomposition F = L*L' (dpptrf)
+   ! (5) compute the reciprocal condition number via dppcon
+   ! (6) Compute P_{t|t} = (I-KZ)P(I-KZ)' + KHK', which ensures symmetry
+
+end subroutine mexFunction

tests/kalman/likelihood/compare_kalman_mex.m

@@ -0,0 +1,116 @@
+function [flag] = compare_kalman_mex(experience)
+
+pp = experience.Number0fObservedVariables;
+mm = experience.SizeOfTheStateVector;
+rr = experience.NumberOfStructuralShocks;
+measurement_error_flag = experience.MeasurementErrors;
+gend = experience.NumberOfPeriods;
+kalman_filter_matlab = experience.kalman_filter_matlab;
+kalman_filter_mex = experience.kalman_filter_mex;
+
+%% SET VARIOUS PARAMETERS 
+kalman_tol = 1e-12;
+riccati_tol =1e-9;
+
+%% SET THE STATE SPACE MODEL:
+
+% I randomly choose the mm eigenvalues of the transition matrix.
+TransitionEigenvalues  = rand(mm,1)*2-1;
+% I randomly choose the mm eigenvectors of the transition matrix
+tmp = rand(mm,mm*100);
+TransitionEigenvectors = tmp*tmp'/(mm*100);
+TransitionEigenvectors = rand(mm,mm);
+% I build the transition matrix
+T = TransitionEigenvectors*diag(TransitionEigenvalues)/TransitionEigenvectors;
+% I randomly choose matrix R
+R = randn(mm,rr);
+% I randomly choose the covariance matrix of the structural innovations
+E = randn(rr,20*rr);
+Q = E*transpose(E)/(20*rr);
+% If needed I randomly choose the covariance matrix of the measurement errors 
+if measurement_error_flag == 0
+    H = zeros(pp,1);
+elseif measurement_error_flag == 1
+    H = rand(pp,1);
+elseif measurement_error_flag == 2
+    E = randn(pp,20*pp);
+    H = E*transpose(E)/(20*pp);
+else
+    disp('compare_kalman_mex: unknown option!')
+end
+% Set the selection vector (mf) 
+MF = transpose(randperm(mm));
+mf = MF(1:pp);
+
+P = lyapunov_symm(T,R*Q*R',riccati_tol,1.000001, riccati_tol);
+
+%% BUILD DATA SET (zero mean):
+a = zeros(mm,1);
+if measurement_error_flag == 0
+    Y = simul_state_space_model(T,R,Q,mf,gend);
+elseif measurement_error_flag == 1
+    H = rand(pp,1);
+    Y = simul_state_space_model(T,R,Q,mf,gend,diag(H));
+elseif measurement_error_flag == 2
+    E = randn(pp,20*pp);
+    H = E*transpose(E)/(20*pp);
+    Y = simul_state_space_model(T,R,Q,mf,gend,H);
+else
+    disp('compare_kalman_mex: unknown option!');
+end
+
+if measurement_error_flag==0
+   HH = 0;
+elseif measurement_error_flag==1
+   HH = diag(H);
+elseif measurement_error_flag==2
+   HH = H;
+end
+
+flag = 0;
+Zflag = 0;
+tic;
+[LIK_matlab,lik_matlab] = kalman_filter_matlab(Y,1,gend,a,P,kalman_tol,riccati_tol,0,0,T,Q,R,HH,mf,mm,pp,rr,Zflag,0,0);
+T_matlab = toc;
+
+tic;
+[LIK_mex,lik_mex] = kalman_filter_mex(Y,a,P,kalman_tol,riccati_tol,T,Q,R,mf,Zflag,HH);
+T_mex = toc;
+
+if T_matlab<T_mex
+   dprintf('Zflag = 0: Matlab Kalman filter is %5.2f times faster than its Fortran counterpart.', T_mex/T_matlab)
+else
+   dprintf('Zflag = 0: Fortran Kalman filter is %5.2f times faster than its Matlab counterpart.', T_matlab/T_mex)
+end
+
+if ((abs((LIK_matlab - LIK_mex)/LIK_matlab) > 1e-6) || (max(abs((lik_matlab-lik_mex)./lik_matlab)) > 1e-6))
+   dprintf("Zflag = 0: discrepancy between Matlab and Fortran Kalman filter results!")
+   flag = 1;
+end
+
+Zflag = 1;
+Z = eye(mm);
+Z = Z(mf,:);
+tic;
+[LIK_matlab_z,lik_matlab_z] = kalman_filter_matlab(Y,1,gend,a,P,kalman_tol,riccati_tol,0,0,T,Q,R,HH,Z,mm,pp,rr,Zflag,0,0);
+T_matlab = toc;
+
+tic;
+[LIK_mex_z,lik_mex_z] = kalman_filter_mex(Y,a,P,kalman_tol,riccati_tol,T,Q,R,Z,Zflag,HH);
+T_mex = toc;
+
+if T_matlab<T_mex
+   dprintf('Zflag = 1: Matlab Kalman filter is %5.2f times faster than its Fortran counterpart.', T_mex/T_matlab)
+else
+   dprintf('Zflag = 1: Fortran Kalman filter is %5.2f times faster than its Matlab counterpart.', T_matlab/T_mex)
+end
+
+if ((abs((LIK_matlab_z - LIK_mex_z)/LIK_matlab_z) > 1e-6) || (max(abs((lik_matlab_z-lik_mex_z)./lik_matlab_z)) > 1e-6))
+   dprintf("Zflag = 1: discrepancy between Matlab and Fortran Kalman filter results!")
+   flag = 1;
+end
+
+if ((abs((LIK_mex - LIK_mex_z)/LIK_mex) > 1e-6) || (max(abs((lik_mex-lik_mex_z)./lik_mex)) > 1e-6))
+   dprintf("Zflag = 1: discrepancy between results with and without the Zflag!")
+   flag = 1;
+end

tests/kalman/likelihood/test_kalman_mex.m

@@ -0,0 +1,108 @@
+debug = true;
+
+if debug
+    source_dir = fileparts(mfilename('fullpath'));
+    source_dir = sprintf('%s/../../../', source_dir)
+else
+    source_dir = getenv('source_root');
+end
+
+mex_path = [source_dir filesep 'build-matlab'];
+addpath(mex_path);
+which kalman_filter
+kalman_filter_mex = @kalman_filter;
+rmpath(mex_path);
+
+matlab_path = [source_dir filesep 'matlab' filesep 'kalman' filesep 'likelihood'];
+addpath(matlab_path);
+which kalman_filter
+kalman_filter_matlab = @kalman_filter;
+rmpath(matlab_path);
+
+addpath([source_dir filesep 'matlab']);
+
+dynare_config;
+
+addpath([source_dir filesep 'tests' filesep 'kalman' filesep 'likelihood']);
+
+testFailed = 0;
+
+if ~debug
+    skipline()
+    disp('***  TESTING: test_kalman_mex.m ***');
+end
+
+t0 = clock;
+
+dprintf('Test 1: No measurement error')
+
+Experience.Number0fObservedVariables = 10;
+Experience.SizeOfTheStateVector = 100;
+Experience.NumberOfStructuralShocks = 12;
+Experience.MeasurementErrors = 0;
+Experience.NumberOfPeriods = 300;
+Experience.kalman_filter_matlab = kalman_filter_matlab;
+Experience.kalman_filter_mex = kalman_filter_mex;
+
+try
+    flag = compare_kalman_mex(Experience);
+    if (flag)
+        testFailed = testFailed+1;
+        if debug
+            dprintf('MEX and MATLAB Kalman filters lead to different results')
+        end
+    end
+catch
+    testFailed = testFailed+1;
+    if debug
+        dprintf('Comparison between MEX and MATLAB Kalman filters failed')
+    end
+end
+
+dprintf('Test 2: measurement error with diagonal variance-covariance matrix')
+Experience.MeasurementErrors = 1;
+
+try
+    flag = compare_kalman_mex(Experience);
+    if (flag)
+        testFailed = testFailed+1;
+        if debug
+            dprintf('MEX and MATLAB Kalman filters lead to different results')
+        end
+    end
+catch
+    testFailed = testFailed+1;
+    if debug
+        dprintf('Comparison between MEX and MATLAB Kalman filters failed')
+    end
+end
+
+dprintf('Test 3: measurement error with general variance-covariance matrix')
+Experience.Number0fObservedVariables = 50;
+Experience.SizeOfTheStateVector = 70;
+Experience.NumberOfStructuralShocks = 50;
+Experience.MeasurementErrors = 2;
+Experience.NumberOfPeriods = 300;
+
+try
+   flag = compare_kalman_mex(Experience);
+   if (flag)
+      testFailed = testFailed+1;
+      if debug
+         dprintf('MEX and MATLAB Kalman filters lead to different results')
+      end
+   end
+catch
+   testFailed = testFailed+1;
+   if debug
+      dprintf('Comparison between MEX and MATLAB Kalman filters failed')
+   end
+end
+
+t1 = clock;
+
+fprintf('\n*** Elapsed time (in seconds): %.1f\n\n', etime(t1, t0));
+
+if ~debug
+    quit(testFailed > 0)
+end