diff --git a/src/filteronestep.f95 b/src/filter1step.f95 similarity index 91% rename from src/filteronestep.f95 rename to src/filter1step.f95 index 8dde156..c98b19e 100644 --- a/src/filteronestep.f95 +++ b/src/filter1step.f95 @@ -1,58 +1,6 @@ -!non-diffuse filtering for single time point -subroutine filteronestep(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,lik,tol,meps,c,p,m,j) - - implicit none - - integer, intent(in) :: p, m,j - integer :: i - integer, intent(in), dimension(p) :: ymiss - double precision, intent(in), dimension(p) :: yt - double precision, intent(in), dimension(m,p) :: zt - double precision, intent(in), dimension(p,p) :: ht - double precision, intent(in), dimension(m,m) :: tt - double precision, dimension(m,m) :: rqr - double precision, intent(in) :: tol,c,meps - double precision, intent(inout) :: lik - double precision, intent(inout), dimension(m) :: at - double precision, intent(inout), dimension(p) :: vt,ft - double precision, intent(inout), dimension(m,p) :: kt - double precision, intent(inout), dimension(m,m) :: pt - double precision, dimension(m,m) :: mm - double precision, dimension(m) :: ahelp - double precision :: finv - - double precision, external :: ddot - - external dgemm, dsymm, dgemv, dsymv, dsyr - - do i = j+1, p - call dsymv('u',m,1.0d0,pt,m,zt(:,i),1,0.0d0,kt(:,i),1) - ft(i) = ddot(m,zt(:,i),1,kt(:,i),1) + ht(i,i) - if(ymiss(i).EQ.0) then - vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1) - if (ft(i) .GT. tol*maxval(zt(:,i)**2)) then - finv = 1.0d0/ft(i) - at = at + vt(i)*finv*kt(:,i) - call dsyr('u',m,-finv,kt(:,i),1,pt,m) - lik = lik - c - 0.5d0*(log(ft(i)) + vt(i)**2*finv) - else - ft(i)=0.0d0 - end if - end if - end do - - call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1) - at = ahelp - - call dsymm('r','u',m,m,1.0d0,pt,m,tt,m,0.0d0,mm,m) - call dgemm('n','t',m,m,m,1.0d0,mm,m,tt,m,0.0d0,pt,m) - pt = pt + rqr - -end subroutine filteronestep - !diffuse filtering for single time point -subroutine diffusefilteronestep(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,& +subroutine dfilter1step(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,& pinf,finf,kinf,rankp,lik,tol,meps,c,p,m,i) implicit none @@ -127,4 +75,57 @@ subroutine diffusefilteronestep(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,& end if end do -end subroutine diffusefilteronestep +end subroutine dfilter1step + +!non-diffuse filtering for single time point +subroutine filter1step(ymiss, yt, zt, ht, tt, rqr, at, pt, vt, ft,kt,lik,tol,c,p,m,j) + + implicit none + + integer, intent(in) :: p, m,j + integer :: i + integer, intent(in), dimension(p) :: ymiss + double precision, intent(in), dimension(p) :: yt + double precision, intent(in), dimension(m,p) :: zt + double precision, intent(in), dimension(p,p) :: ht + double precision, intent(in), dimension(m,m) :: tt + double precision, dimension(m,m) :: rqr + double precision, intent(in) :: tol,c + double precision, intent(inout) :: lik + double precision, intent(inout), dimension(m) :: at + double precision, intent(inout), dimension(p) :: vt,ft + double precision, intent(inout), dimension(m,p) :: kt + double precision, intent(inout), dimension(m,m) :: pt + double precision, dimension(m,m) :: mm + double precision, dimension(m) :: ahelp + double precision :: finv + + double precision, external :: ddot + + external dgemm, dsymm, dgemv, dsymv, dsyr + + do i = j+1, p + call dsymv('u',m,1.0d0,pt,m,zt(:,i),1,0.0d0,kt(:,i),1) + ft(i) = ddot(m,zt(:,i),1,kt(:,i),1) + ht(i,i) + if(ymiss(i).EQ.0) then + vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1) + if (ft(i) .GT. tol*maxval(zt(:,i)**2)) then + finv = 1.0d0/ft(i) + at = at + vt(i)*finv*kt(:,i) + call dsyr('u',m,-finv,kt(:,i),1,pt,m) + lik = lik - c - 0.5d0*(log(ft(i)) + vt(i)**2*finv) + else + ft(i)=0.0d0 + end if + end if + end do + + call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1) + at = ahelp + + call dsymm('r','u',m,m,1.0d0,pt,m,tt,m,0.0d0,mm,m) + call dgemm('n','t',m,m,m,1.0d0,mm,m,tt,m,0.0d0,pt,m) + pt = pt + rqr + +end subroutine filter1step + diff --git a/src/filteronestepnovar.f95 b/src/filter1stepnovar.f95 similarity index 79% rename from src/filteronestepnovar.f95 rename to src/filter1stepnovar.f95 index 7bf86a4..9cfccdf 100644 --- a/src/filteronestepnovar.f95 +++ b/src/filter1stepnovar.f95 @@ -1,75 +1,79 @@ -!non-diffuse filtering for single time point -subroutine filteronestepnv(ymiss, yt, zt, tt, at,vt, ft,kt,p,m,j) +!diffuse filtering for single time point +subroutine dfilter1stepnv(ymiss, yt, zt, tt, at, vt, ft,kt,& +finf,kinf,p,m,j,lik) implicit none - integer, intent(in) :: p, m,j - integer :: i + integer, intent(in) :: p, m, j + integer :: i integer, intent(in), dimension(p) :: ymiss double precision, intent(in), dimension(p) :: yt double precision, intent(in), dimension(m,p) :: zt double precision, intent(in), dimension(m,m) :: tt double precision, intent(inout), dimension(m) :: at - double precision, intent(in), dimension(p) :: ft double precision, intent(inout), dimension(p) :: vt - double precision, intent(in), dimension(m,p) :: kt + double precision, intent(in), dimension(p) :: ft,finf + double precision, intent(in), dimension(m,p) :: kt,kinf + double precision, intent(inout) :: lik double precision, dimension(m) :: ahelp - double precision, external :: ddot - + double precision, external :: ddot external dgemv - do i = j+1, p + do i = 1, j if(ymiss(i).EQ.0) then vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1) - if (ft(i) .GT. 0.0d0) then - at = at + vt(i)/ft(i)*kt(:,i) + if (finf(i) .GT. 0.0d0) then + at = at + vt(i)/finf(i)*kinf(:,i) + lik = lik - 0.5d0*log(finf(i)) + else + if (ft(i) .GT. 0.0d0) then + at = at + vt(i)/ft(i)*kt(:,i) + lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i)) + end if end if end if end do - - call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1) - at = ahelp - -end subroutine filteronestepnv + if(j .EQ. p) then + call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1) + at = ahelp + end if +end subroutine dfilter1stepnv -!diffuse filtering for single time point -subroutine diffusefilteronestepnv(ymiss, yt, zt, tt, at, vt, ft,kt,& -finf,kinf,p,m,j) +!non-diffuse filtering for single time point +subroutine filter1stepnv(ymiss, yt, zt, tt, at,vt, ft,kt,p,m,j,lik) implicit none - integer, intent(in) :: p, m, j - integer :: i + integer, intent(in) :: p, m,j + integer :: i integer, intent(in), dimension(p) :: ymiss double precision, intent(in), dimension(p) :: yt double precision, intent(in), dimension(m,p) :: zt double precision, intent(in), dimension(m,m) :: tt double precision, intent(inout), dimension(m) :: at + double precision, intent(in), dimension(p) :: ft double precision, intent(inout), dimension(p) :: vt - double precision, intent(in), dimension(p) :: ft,finf - double precision, intent(in), dimension(m,p) :: kt,kinf + double precision, intent(in), dimension(m,p) :: kt + double precision, intent(inout) :: lik double precision, dimension(m) :: ahelp double precision, external :: ddot external dgemv - do i = 1, j + do i = j+1, p if(ymiss(i).EQ.0) then vt(i) = yt(i) - ddot(m,zt(:,i),1,at,1) - if (finf(i) .GT. 0.0d0) then - at = at + vt(i)/finf(i)*kinf(:,i) - else - if (ft(i) .GT. 0.0d0) then - at = at + vt(i)/ft(i)*kt(:,i) - end if + if (ft(i) .GT. 0.0d0) then + at = at + vt(i)/ft(i)*kt(:,i) + lik = lik - 0.5d0*(log(ft(i)) + vt(i)**2/ft(i)) end if end if end do - if(j .EQ. p) then - call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1) - at = ahelp - end if -end subroutine diffusefilteronestepnv + + call dgemv('n',m,m,1.0d0,tt,m,at,1,0.0d0,ahelp,1) + at = ahelp + +end subroutine filter1stepnv diff --git a/src/filtersimfast.f95 b/src/filtersimfast.f95 index 9f27766..64b2511 100644 --- a/src/filtersimfast.f95 +++ b/src/filtersimfast.f95 @@ -1,11 +1,11 @@ !fast filtering algorithm used in simulation filter subroutine filtersimfast(yt, ymiss, timevar, zt,tt, & -a1, ft,kt,finf, kinf, dt, jt, p, m, n,tol,at) +a1, ft,kt,finf, kinf, dt, jt, p, m, n,at) implicit none integer, intent(in) :: p, m,n,dt,jt - integer :: t, j + integer :: t integer, intent(in), dimension(n,p) :: ymiss integer, intent(in), dimension(5) :: timevar double precision, intent(in), dimension(n,p) :: yt @@ -14,39 +14,42 @@ subroutine filtersimfast(yt, ymiss, timevar, zt,tt, & double precision, intent(in), dimension(m) :: a1 double precision, intent(in), dimension(p,n) :: ft,finf double precision, intent(in), dimension(m,p,n) :: kt,kinf - double precision, intent(in) :: tol double precision, intent(inout), dimension(m,n+1) :: at double precision, dimension(p,n) :: vt - double precision, dimension(m) :: arec + double precision :: lik double precision, external :: ddot external dgemv + lik = 0.0d0 + at(:,1) = a1 - !diffuse filtering begins - do t = 1, (dt - 1) + if(dt .GT. 0) then + !diffuse filtering begins + do t = 1, (dt - 1) + at(:,t+1) = at(:,t) + call dfilter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),& + finf(:,t),kinf(:,:,t),p,m,p,lik) + end do + + t = dt at(:,t+1) = at(:,t) - call diffusefilteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + call dfilter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),& - finf(:,t),kinf(:,:,t),p,m,p) + finf(:,t),kinf(:,:,t),p,m,jt,lik) + !non-diffuse filtering begins + if(jt .LT. p) then + call filter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),p,m,jt,lik) + end if + end if + !Non-diffuse filtering continues from t=d+1, i=1 + do t = dt + 1, n + at(:,t+1) = at(:,t) + call filter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),p,m,0,lik) end do - t = dt - at(:,t+1) = at(:,t) - call diffusefilteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& - tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),& - finf(:,t),kinf(:,:,t),p,m,jt) - !non-diffuse filtering begins - call filteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& - tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),p,m,jt) - - if(dt.LT.n) then - !Non-diffuse filtering continues from t=d+1, i=1 - do t = dt + 1, n - at(:,t+1) = at(:,t) - call filteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& - tt(:,:,(t-1)*timevar(3)+1),at(:,t+1),vt(:,t),ft(:,t),kt(:,:,t),p,m,0) - end do - end if end subroutine filtersimfast diff --git a/src/gloglik.f95 b/src/gloglik.f95 index f0f91d4..2161b0b 100644 --- a/src/gloglik.f95 +++ b/src/gloglik.f95 @@ -8,7 +8,7 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,& integer, intent(in) :: p, m, r, n integer, intent(inout) :: rankp,marginal - integer :: t, i,d,j,tv + integer :: t,d,j,tv integer, intent(in), dimension(n,p) :: ymiss integer, intent(in), dimension(5) :: timevar double precision, intent(in), dimension(n,p) :: yt @@ -21,12 +21,12 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,& double precision, intent(in), dimension(m,m) :: p1,p1inf double precision, intent(in) :: tol double precision, intent(inout) :: lik - double precision, dimension(m) :: at,arec + double precision, dimension(m) :: at double precision, dimension(p) :: vt,ft,finf double precision, dimension(m,p) :: kt,kinf - double precision, dimension(m,m) :: pt,pinf,mm + double precision, dimension(m,m) :: pt,pinf double precision, dimension(m,r) :: mr - double precision :: c, meps, finv + double precision :: c, meps double precision, external :: ddot double precision, dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr @@ -56,15 +56,15 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,& diffuse: do while(d .LT. n .AND. rankp .GT. 0) d = d+1 - call diffusefilteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& + call dfilter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1), at,pt,vt,ft,kt,pinf,finf,kinf,rankp,lik,tol,meps,c,p,m,j) end do diffuse if(rankp .EQ. 0 .AND. j .LT. p) then !non-diffuse filtering begins - call filteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& - tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),at,pt,vt,ft,kt,lik,tol,meps,c,p,m,j) + call filter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& + tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),at,pt,vt,ft,kt,lik,tol,c,p,m,j) else j = p @@ -77,8 +77,8 @@ subroutine gloglik(yt, ymiss, timevar, zt, ht, tt, rt, qt, a1, p1, p1inf,& !Non-diffuse filtering continues from t=d+1, i=1 do t = d+1, n - call filteronestep(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),& - tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),at,pt,vt,ft,kt,lik,tol,meps,c,p,m,0) + call filter1step(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),& + tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),at,pt,vt,ft,kt,lik,tol,c,p,m,0) end do diff --git a/src/kfilter.f95 b/src/kfilter.f95 index 7147f89..5bb5c4d 100644 --- a/src/kfilter.f95 +++ b/src/kfilter.f95 @@ -7,7 +7,7 @@ subroutine kfilter(yt, ymiss, timevar, zt, ht,tt, rt, qt, a1, p1, p1inf, p,n,m,r integer, intent(in) :: p, m, r, n,filtersignal integer, intent(inout) :: d, j, rankp - integer :: t, i,tv + integer :: t,tv integer, intent(in), dimension(n,p) :: ymiss integer, intent(in), dimension(5) :: timevar double precision, intent(in), dimension(n,p) :: yt @@ -26,11 +26,9 @@ subroutine kfilter(yt, ymiss, timevar, zt, ht,tt, rt, qt, a1, p1, p1inf, p,n,m,r double precision, intent(inout) :: lik double precision, intent(inout), dimension(p,p,n) :: thetavar double precision, intent(inout), dimension(n,p) :: theta - double precision, dimension(m) :: arec - double precision, dimension(m,m) ::pirec,mm,prec double precision, dimension(m,r) :: mr double precision, dimension(p,m) :: pm - double precision :: c,meps,finv + double precision :: c,meps double precision, external :: ddot double precision, dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr external dgemm, dsymm, dgemv, dsymv, dsyr, dsyr2 @@ -49,10 +47,7 @@ subroutine kfilter(yt, ymiss, timevar, zt, ht,tt, rt, qt, a1, p1, p1inf, p,n,m,r d=0 pinf(:,:,1)=p1inf pt(:,:,1) = p1 - prec = p1 - pirec = p1inf at(:,1) = a1 - arec = a1 ! diffuse initialization if(rankp .GT. 0) then diffuse: do while(d .LT. n .AND. rankp .GT. 0) @@ -60,7 +55,7 @@ subroutine kfilter(yt, ymiss, timevar, zt, ht,tt, rt, qt, a1, p1, p1inf, p,n,m,r at(:,d+1) = at(:,d) pt(:,:,d+1) = pt(:,:,d) pinf(:,:,d+1) = pinf(:,:,d) - call diffusefilteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& + call dfilter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),& at(:,d+1),pt(:,:,d+1),vt(:,d),ft(:,d),kt(:,:,d),pinf(:,:,d+1),finf(:,d),kinf(:,:,d),rankp,lik,tol,meps,c,p,m,j) end do diffuse @@ -69,9 +64,9 @@ subroutine kfilter(yt, ymiss, timevar, zt, ht,tt, rt, qt, a1, p1, p1inf, p,n,m,r if(rankp .EQ. 0 .AND. j .LT. p) then !non-diffuse filtering begins - call filteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& + call filter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),& - at(:,d+1),pt(:,:,d+1),vt(:,d),ft(:,d),kt(:,:,d),lik,tol,meps,c,p,m,j) + at(:,d+1),pt(:,:,d+1),vt(:,d),ft(:,d),kt(:,:,d),lik,tol,c,p,m,j) else j = p @@ -84,9 +79,9 @@ subroutine kfilter(yt, ymiss, timevar, zt, ht,tt, rt, qt, a1, p1, p1inf, p,n,m,r do t = d+1, n at(:,t+1) = at(:,t) pt(:,:,t+1) = pt(:,:,t) - call filteronestep(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),& + call filter1step(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),& tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),& - at(:,t+1),pt(:,:,t+1),vt(:,t),ft(:,t),kt(:,:,t),lik,tol,meps,c,p,m,0) + at(:,t+1),pt(:,:,t+1),vt(:,t),ft(:,t),kt(:,:,t),lik,tol,c,p,m,0) end do if(filtersignal.EQ.1) then diff --git a/src/kfstheta.f95 b/src/kfstheta.f95 index a4f6b25..d4ffc89 100644 --- a/src/kfstheta.f95 +++ b/src/kfstheta.f95 @@ -21,43 +21,43 @@ subroutine kfstheta(yt, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, tv, a1, p1, p1inf double precision, dimension(p,n) :: ft,finf double precision, dimension(m,p,n) :: kt,kinf double precision, dimension(p,n) :: vt - double precision, dimension(m) :: at, arec - double precision, dimension(m,m) :: mm,pinf,pt + double precision, dimension(m) :: at,mhelp + double precision, dimension(m,m) :: pinf,pt double precision, intent(in) :: tol - double precision, dimension(m) :: rrec,rrec1,rhelp,help - double precision, dimension(m,m) :: im,linf,l0 + double precision, dimension(m) :: rt0,rt1 + double precision, dimension(m,m) :: im double precision, dimension(r,n) :: etahat - double precision :: c, meps, finv + double precision :: c, meps double precision, external :: ddot double precision, intent(inout), dimension(n,p) :: thetahat double precision, intent(inout) :: lik - external dgemm, dsymm, dgemv, dsymv, dsyr, dsyr2, dger + double precision, dimension(p) :: epshat + external dgemv, dsymv + + epshat = 0.0d0 meps = epsilon(meps) lik=0.0d0 c = 0.5d0*log(8.0d0*atan(1.0d0)) - im = 0.0d0 - do i = 1, m - im(i,i) = 1.0d0 - end do j=0 d=0 pt = p1 pinf = p1inf at = a1 if(rankp .GT. 0) then + !diffuse filtering diffuse: do while(d .LT. n .AND. rankp .GT. 0) d = d+1 - call diffusefilteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& + call dfilter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),& at,pt,vt(:,d),ft(:,d),kt(:,:,d),pinf,finf(:,d),kinf(:,:,d),rankp,lik,tol,meps,c,p,m,j) end do diffuse if(rankp .EQ. 0 .AND. j .LT. p) then - call filteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& + call filter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),& - at,pt,vt(:,d),ft(:,d),kt(:,:,d),lik,tol,meps,c,p,m,j) + at,pt,vt(:,d),ft(:,d),kt(:,:,d),lik,tol,c,p,m,j) else j = p @@ -67,147 +67,55 @@ subroutine kfstheta(yt, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, tv, a1, p1, p1inf !Non-diffuse filtering continues from t=d+1, i=1 do t = d+1, n - call filteronestep(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),& + call filter1step(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),& tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),& - at,pt,vt(:,t),ft(:,t),kt(:,:,t),lik,tol,meps,c,p,m,0) + at,pt,vt(:,t),ft(:,t),kt(:,:,t),lik,tol,c,p,m,0) end do - !smoothing begins + !smoothing begins + im = 0.0d0 + do i = 1, m + im(i,i) = 1.0d0 + end do - rrec = 0.0d0 + rt0 = 0.0d0 do t = n, d+1, -1 - call dgemv('t',m,r,1.0d0,rtv(:,:,(t-1)*timevar(4)+1),m,rrec,1,0.0d0,help,1) - call dsymv('l',r,1.0d0,qt(:,:,(t-1)*timevar(5)+1),r,help,1,0.0d0,etahat(:,t),1) - call dgemv('t',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,rrec,1,0.0d0,rhelp,1) - rrec = rhelp - do i = p, 1 , -1 - if(ymiss(t,i).EQ.0) then - if(ft(i,t) .GT. 0.0d0) then - finv = 1.0d0/ft(i,t) - l0 = im - call dger(m,m,-finv,kt(:,i,t),1,zt(i,:,(t-1)*timevar(1)+1),1,l0,m) - call dgemv('t',m,m,1.0d0,l0,m,rrec,1,0.0d0,rhelp,1) - rrec = rhelp + vt(i,t)*finv*zt(i,:,(t-1)*timevar(1)+1) - end if - end if - end do + call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & + ft(:,t),kt(:,:,t), im,p,m,r,1,rt0,etahat(:,t),epshat,.FALSE.) end do - if(d.GT.0) then - t=d - call dgemv('t',m,r,1.0d0,rtv(:,:,(t-1)*timevar(4)+1),m,rrec,1,0.0d0,help,1) - call dsymv('l',r,1.0d0,qt(:,:,(t-1)*timevar(5)+1),r,help,1,0.0d0,etahat(:,t),1) - call dgemv('t',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,rrec,1,0.0d0,rhelp,1) - rrec = rhelp - do i = p, (j+1) , -1 - if(ymiss(t,i).EQ.0) then - if(ft(i,t) .GT. 0.0d0) then - finv = 1.0d0/ft(i,t) - l0 = im - call dger(m,m,-finv,kt(:,i,t),1,zt(i,:,(t-1)*timevar(1)+1),1,l0,m) - call dgemv('t',m,m,1.0d0,l0,m,rrec,1,0.0d0,rhelp,1) - rrec = rhelp + vt(i,t)*finv*zt(i,:,(t-1)*timevar(1)+1) - end if - end if - end do - rrec1 = 0.0d0 - do i = j, 1, -1 - if(ymiss(t,i).EQ.0) then - if(finf(i,t).GT. 0.0d0) then - finv = 1.0d0/finf(i,t) - - linf = im - call dger(m,m,-finv,kinf(:,i,t),1,zt(i,:,(t-1)*timevar(1)+1),1,linf,m) - - rhelp = kinf(:,i,t)*ft(i,t)*finv - kt(:,i,t) - l0=0.0d0 - call dger(m,m,finv,rhelp,1,zt(i,:,(t-1)*timevar(1)+1),1,l0,m) - - call dgemv('t',m,m,1.0d0,linf,m,rrec1,1,0.0d0,rhelp,1) - rrec1 = rhelp - call dgemv('t',m,m,1.0d0,l0,m,rrec,1,1.0d0,rrec1,1) - rrec1 = rrec1 + vt(i,t)*finv*zt(i,:,(t-1)*timevar(1)+1) - - call dgemv('t',m,m,1.0d0,linf,m,rrec,1,0.0d0,rhelp,1) - rrec = rhelp - else - if(ft(i,t).GT. 0.0d0) then - finv = 1.0d0/ft(i,t) - l0 = im - call dger(m,m,-finv,kt(:,i,t),1,zt(i,:,(t-1)*timevar(1)+1),1,l0,m) - - call dgemv('t',m,m,1.0d0,l0,m,rrec,1,0.0d0,rhelp,1) - rrec = rhelp + vt(i,t)*finv*zt(i,:,(t-1)*timevar(1)+1) - - call dgemv('t',m,m,1.0d0,l0,m,rrec1,1,0.0d0,rhelp,1) - rrec1 = rhelp - end if - end if - end if - end do - - do t=(d-1), 1, -1 - call dgemv('t',m,r,1.0d0,rtv(:,:,(t-1)*timevar(4)+1),m,rrec,1,0.0d0,help,1) - call dsymv('l',r,1.0d0,qt(:,:,(t-1)*timevar(5)+1),r,help,1,0.0d0,etahat(:,t),1) - call dgemv('t',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,rrec,1,0.0d0,rhelp,1) - rrec = rhelp - call dgemv('t',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,rrec1,1,0.0d0,rhelp,1) - rrec1 = rhelp - - do i = p, 1, -1 - if(ymiss(t,i).EQ.0) then - if(finf(i,t).GT. 0.0d0) then - finv = 1.0d0/finf(i,t) - - linf = im - call dger(m,m,-finv,kinf(:,i,t),1,zt(i,:,(t-1)*timevar(1)+1),1,linf,m) - - rhelp = kinf(:,i,t)*ft(i,t)*finv - kt(:,i,t) - l0=0.0d0 - call dger(m,m,finv,rhelp,1,zt(i,:,(t-1)*timevar(1)+1),1,l0,m) - - call dgemv('t',m,m,1.0d0,linf,m,rrec1,1,0.0d0,rhelp,1) - rrec1 = rhelp - call dgemv('t',m,m,1.0d0,l0,m,rrec,1,1.0d0,rrec1,1) - rrec1 = rrec1 + vt(i,t)*finv*zt(i,:,(t-1)*timevar(1)+1) - - call dgemv('t',m,m,1.0d0,linf,m,rrec,1,0.0d0,rhelp,1) - rrec = rhelp - else - if(ft(i,t).GT. 0.0d0) then - finv = 1.0d0/ft(i,t) - l0 = im - call dger(m,m,-finv,kt(:,i,t),1,zt(i,:,(t-1)*timevar(1)+1),1,l0,m) - - call dgemv('t',m,m,1.0d0,l0,m,rrec,1,0.0d0,rhelp,1) - rrec = rhelp + vt(i,t)*finv*zt(i,:,(t-1)*timevar(1)+1) - - call dgemv('t',m,m,1.0d0,l0,m,rrec1,1,0.0d0,rhelp,1) - rrec1 = rhelp - - end if - - end if - end if - end do - + if(d .GT. 0) then + t = d + if(j .LT. p) then + call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & + ft(:,t),kt(:,:,t), im,p,m,r,j+1,rt0,etahat(:,t),epshat,.FALSE.) + end if + rt1 = 0.0d0 + call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & + ft(:,t),kt(:,:,t), im,p,m,r,j,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat,.FALSE.) + do t = (d - 1), 1, -1 + call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & + ft(:,t),kt(:,:,t), im,p,m,r,p,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat,.FALSE.) end do end if - at = a1 - call dsymv('l',m,1.0d0,p1,m,rrec,1,1.0d0,at,1) + call dsymv('l',m,1.0d0,p1,m,rt0,1,1.0d0,at,1) if(d .GT. 0) then - call dsymv('l',m,1.0d0,p1inf,m,rrec1,1,1.0d0,at,1) + call dsymv('l',m,1.0d0,p1inf,m,rt1,1,1.0d0,at,1) end if call dgemv('n',p,m,1.0d0,zt(:,:,1),p,at,1,0.0d0,thetahat(1,:),1) do t = 2, n - call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,at,1,0.0d0,help,1) - at=help + call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,at,1,0.0d0,mhelp,1) + at=mhelp call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etahat(:,t-1),1,1.0d0,at,1) call dgemv('n',p,m,1.0d0,zt(:,:,(t-1)*timevar(1)+1),p,at,1,0.0d0,thetahat(t,:),1) diff --git a/src/ptheta.f95 b/src/ptheta.f95 index 38560f9..ff0760e 100644 --- a/src/ptheta.f95 +++ b/src/ptheta.f95 @@ -1,6 +1,5 @@ ! functions for computing p(theta) - ! Only differences with gloglik is what is stored/returned (ft,finf,kt,kinf) ! Also no missing values as yt is actually the signals, and no ht subroutine pthetafirst(yt, timevar, zt, tt, rqr, a1, p1, p1inf,& @@ -11,7 +10,7 @@ subroutine pthetafirst(yt, timevar, zt, tt, rqr, a1, p1, p1inf,& integer, intent(in) :: p, m, n integer, intent(inout) :: rankp2,d,j - integer :: t, i,tv,rankp + integer :: t, tv,rankp integer, intent(in), dimension(5) :: timevar double precision, intent(in), dimension(n,p) :: yt double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt @@ -20,14 +19,13 @@ subroutine pthetafirst(yt, timevar, zt, tt, rqr, a1, p1, p1inf,& double precision, intent(in), dimension(m,m) :: p1,p1inf double precision, intent(in) :: tol double precision, intent(inout) :: lik - double precision, dimension(m) :: at,arec + double precision, dimension(m) :: at double precision, dimension(p) :: vt double precision, intent(inout), dimension(p,n) :: ft,finf double precision, intent(inout), dimension(m,p,n) :: kt,kinf - double precision, dimension(m,m) :: pt,pinf,mm - double precision, external :: ddot + double precision, dimension(m,m) :: pt,pinf double precision, intent(inout), dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr - double precision :: meps, finv + double precision :: meps integer, dimension(p) :: ymiss double precision, dimension(p,p) :: ht @@ -49,31 +47,25 @@ subroutine pthetafirst(yt, timevar, zt, tt, rqr, a1, p1, p1inf,& if(rankp .GT. 0) then diffuse: do while(d .LT. n .AND. rankp .GT. 0) d = d+1 - call diffusefilteronestep(ymiss,yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht,& + call dfilter1step(ymiss,yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht,& tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),& at,pt,vt,ft(:,d),kt(:,:,d),pinf,finf(:,d),kinf(:,:,d),rankp,lik,tol,meps,0.0d0,p,m,j) - - end do diffuse if(rankp .EQ. 0 .AND. j .LT. p) then - call filteronestep(ymiss,yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht,& + call filter1step(ymiss,yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht,& tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),& - at,pt,vt,ft(:,d),kt(:,:,d),lik,tol,meps,0.0d0,p,m,j) - + at,pt,vt,ft(:,d),kt(:,:,d),lik,tol,0.0d0,p,m,j) else j = p end if - end if - - !Non-diffuse filtering continues from t=d+1, i=1 do t = d+1, n - call filteronestep(ymiss,yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht,& + call filter1step(ymiss,yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht,& tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),& - at,pt,vt,ft(:,t),kt(:,:,t),lik,tol,meps,0.0d0,p,m,0) + at,pt,vt,ft(:,t),kt(:,:,t),lik,tol,0.0d0,p,m,0) end do @@ -90,106 +82,46 @@ subroutine pthetarest(yt, timevar, zt, tt, a1,& implicit none integer, intent(in) :: p, m, n,dt,jt - integer :: i,j,t,d + integer :: t integer, intent(in), dimension(5) :: timevar double precision, intent(in), dimension(n,p) :: yt double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt double precision, intent(in), dimension(m) :: a1 double precision, intent(inout) :: lik - double precision, dimension(m) :: at,arec + double precision, dimension(m) :: at double precision, dimension(p) :: vt double precision, intent(in), dimension(p,n) :: ft,finf double precision, intent(in), dimension(m,p,n) :: kt,kinf + integer, dimension(p) :: ymiss double precision, external :: ddot external dgemv - - j=0 - d=0 - arec = a1 - - if(dt.GT.0) then - diffuse: do while(d .LT. (dt-1)) - d = d+1 - do j=1, p - - vt(j) = yt(d,j) - ddot(m,zt(j,:,(d-1)*timevar(1)+1),1,arec,1) - if (finf(j,d) .GT. 0.0d0) then - arec = arec + vt(j)/finf(j,d)*kinf(:,j,d) - - lik = lik - 0.5d0*log(finf(j,d)) - else - if(ft(j,d) .GT. 0.0d0) then - arec = arec + vt(j)/ft(j,d)*kt(:,j,d) - - lik = lik - 0.5d0*(log(ft(j,d)) + vt(j)**2/ft(j,d)) - end if - end if - - end do - - call dgemv('n',m,m,1.0d0,tt(:,:,(d-1)*timevar(3)+1),m,arec,1,0.0d0,at,1) - arec = at - - end do diffuse - - d = dt - do j=1, jt - - vt(j) = yt(d,j) - ddot(m,zt(j,:,(d-1)*timevar(1)+1),1,arec,1) - - if (finf(j,d) .GT. 0.0d0) then - arec = arec + vt(j)/finf(j,d)*kinf(:,j,d) - lik = lik - 0.5d0*log(finf(j,d)) - else - if(ft(j,d) .GT. 0.0d0) then - arec = arec + vt(j)/ft(j,d)*kt(:,j,d) - lik = lik - 0.5d0*(log(ft(j,d)) + vt(j)**2/ft(j,d)) - end if - end if - + ymiss = 0 + at = a1 + if(dt .GT. 0) then + !diffuse filtering begins + do t = 1, dt - 1 + call dfilter1stepnv(ymiss,yt(t,:),& + transpose(zt(:,:,(t-1)*timevar(1)+1)),tt(:,:,(t-1)*timevar(3)+1),& + at,vt,ft(:,t),kt(:,:,t), finf(:,t),kinf(:,:,t),p,m,p,lik) end do - + t = dt + call dfilter1stepnv(ymiss,yt(t,:),& + transpose(zt(:,:,(t-1)*timevar(1)+1)),tt(:,:,(t-1)*timevar(3)+1),& + at,vt,ft(:,t),kt(:,:,t),finf(:,t),kinf(:,:,t),p,m,jt,lik) !non-diffuse filtering begins - - do i = jt+1, p - - vt(i) = yt(d,i) - ddot(m,zt(i,:,(d-1)*timevar(1)+1),1,arec,1) - if (ft(i,d) .GT. 0.0d0) then !ft.NE.0 - arec = arec + vt(i)/ft(i,d)*kt(:,i,d) - - lik = lik - 0.5d0*(log(ft(i,d)) + vt(i)**2/ft(i,d)) - end if - - end do - - call dgemv('n',m,m,1.0d0,tt(:,:,(d-1)*timevar(3)+1),m,arec,1,0.0d0,at,1) - arec = at - end if - - if(dt.LT.n) then - - !Non-diffuse filtering continues from t=d+1, i=1 - do t = dt+1, n - do i = 1, p - - vt(i) = yt(t,i) - ddot(m,zt(i,:,(t-1)*timevar(1)+1),1,arec,1) - if (ft(i,t) .GT. 0.0d0) then - arec = arec + vt(i)/ft(i,t)*kt(:,i,t) - lik = lik - 0.5d0*(log(ft(i,t)) + vt(i)**2/ft(i,t)) - end if - - end do - - call dgemv('n',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,arec,1,0.0d0,at,1) - arec = at - end do - + if(jt .LT. p) then + call filter1stepnv(ymiss,yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + tt(:,:,(t-1)*timevar(3)+1),at,vt,ft(:,t),kt(:,:,t),p,m,jt,lik) + end if end if - - + !Non-diffuse filtering continues from t=d+1, i=1 + do t = dt + 1, n + call filter1stepnv(ymiss,yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + tt(:,:,(t-1)*timevar(3)+1),at,vt,ft(:,t),kt(:,:,t),p,m,0,lik) + end do end subroutine pthetarest diff --git a/src/simfilter.f95 b/src/simfilter.f95 index 4dcae6d..ea78cc6 100644 --- a/src/simfilter.f95 +++ b/src/simfilter.f95 @@ -116,7 +116,7 @@ subroutine simfilter(ymiss,timevar, yt, zt, ht, tt, rtv, qt, a1, p1, & atplus=0.0d0 call filtersimfast(yplus, ymiss, timevar, zt,tt, a1, ft,kt,& - finf, kinf, d, j, p, m, n,tol,atplus) + finf, kinf, d, j, p, m, n,atplus) if(simwhat.EQ.4) then do t = 1, n sim(:,t,i) = at(:,t) - atplus(:,t) + aplus(:,t) diff --git a/src/smoothonestep.f95 b/src/smoothonestep.f95 index 1bc8541..8da4682 100644 --- a/src/smoothonestep.f95 +++ b/src/smoothonestep.f95 @@ -1,5 +1,5 @@ !non-diffuse smoothing for single time point -subroutine smoothonestep(ymiss, zt, ht, tt, rtv, qt, vt, ft,kt,& +subroutine smooth1step(ymiss, zt, ht, tt, rtv, qt, vt, ft,kt,& im,p,m,r,j,rt,etahat,epshat,needeps) implicit none @@ -47,9 +47,9 @@ subroutine smoothonestep(ymiss, zt, ht, tt, rtv, qt, vt, ft,kt,& end if end do -end subroutine smoothonestep +end subroutine smooth1step -subroutine diffusesmoothonestep(ymiss, zt, ht, tt, rtv, qt, vt, ft,kt,& +subroutine dsmooth1step(ymiss, zt, ht, tt, rtv, qt, vt, ft,kt,& im,p,m,r,j,rt,rt1,finf,kinf,etahat,epshat,needeps) implicit none @@ -128,4 +128,4 @@ subroutine diffusesmoothonestep(ymiss, zt, ht, tt, rtv, qt, vt, ft,kt,& end if end do -end subroutine diffusesmoothonestep +end subroutine dsmooth1step diff --git a/src/smoothsim.f95 b/src/smoothsim.f95 index ddbce63..c54bb61 100644 --- a/src/smoothsim.f95 +++ b/src/smoothsim.f95 @@ -22,15 +22,14 @@ subroutine smoothsim(yt, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, a1, p1, p1inf, & double precision, intent(inout), dimension(p,n) :: ft,finf double precision, intent(inout), dimension(m,p,n) :: kt,kinf double precision, dimension(p,n) :: vt - double precision, dimension(m) :: at, arec - double precision, dimension(m,m) :: pt,pinf,mm + double precision, dimension(m) :: at + double precision, dimension(m,m) :: pt,pinf double precision, intent(in) :: tol - double precision, dimension(m) :: rrec,rrec1,rhelp,help - double precision, dimension(m,m) :: im,linf,l0 + double precision, dimension(m,m) :: im double precision, intent(inout), dimension(r,n) :: etahat double precision, intent(inout), dimension(p,n) :: epshat double precision, intent(inout), dimension(m) :: rt0,rt1 - double precision :: meps, finv,lik + double precision :: meps,lik double precision, external :: ddot external dgemm, dsymm, dgemv, dsymv, dsyr, dsyr2, dger @@ -38,10 +37,7 @@ subroutine smoothsim(yt, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, a1, p1, p1inf, & meps = epsilon(meps) tv = max(timevar(4),timevar(5)) lik = 0.0d0 - im = 0.0d0 - do i = 1, m - im(i,i) = 1.0d0 - end do + j=0 d=0 pinf = p1inf @@ -50,63 +46,57 @@ subroutine smoothsim(yt, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, a1, p1, p1inf, & if(rankp .GT. 0) then diffuse: do while(d .LT. n .AND. rankp .GT. 0) d = d+1 - call diffusefilteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& + call dfilter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),& at,pt,vt(:,d),ft(:,d),kt(:,:,d),pinf,finf(:,d),kinf(:,:,d),rankp,lik,tol,meps,0.0d0,p,m,j) end do diffuse if(rankp .EQ. 0 .AND. j .LT. p) then - call filteronestep(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& + call filter1step(ymiss(d,:),yt(d,:),transpose(zt(:,:,(d-1)*timevar(1)+1)),ht(:,:,(d-1)*timevar(2)+1),& tt(:,:,(d-1)*timevar(3)+1),rqr(:,:,(d-1)*tv+1),& - at,pt,vt(:,d),ft(:,d),kt(:,:,d),lik,tol,meps,0.0d0,p,m,j) - + at,pt,vt(:,d),ft(:,d),kt(:,:,d),lik,tol,0.0d0,p,m,j) else j = p end if - end if !Non-diffuse filtering continues from t=d+1, i=1 do t = d+1, n - call filteronestep(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),& + call filter1step(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),ht(:,:,(t-1)*timevar(2)+1),& tt(:,:,(t-1)*timevar(3)+1),rqr(:,:,(t-1)*tv+1),& - at,pt,vt(:,t),ft(:,t),kt(:,:,t),lik,tol,meps,0.0d0,p,m,0) + at,pt,vt(:,t),ft(:,t),kt(:,:,t),lik,tol,0.0d0,p,m,0) end do !smoothing begins + im = 0.0d0 + do i = 1, m + im(i,i) = 1.0d0 + end do + rt0 = 0.0d0 do t = n, d+1, -1 - call smoothonestep(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & ft(:,t),kt(:,:,t), im,p,m,r,1,rt0,etahat(:,t),epshat(:,t),needeps) - end do if(d .GT. 0) then t = d if(j .LT. p) then - call smoothonestep(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & ft(:,t),kt(:,:,t), im,p,m,r,j+1,rt0,etahat(:,t),epshat(:,t),needeps) end if - - rt1 = 0.0d0 - call diffusesmoothonestep(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & ft(:,t),kt(:,:,t), im,p,m,r,j,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat(:,t),needeps) - - do t=(d-1), 1, -1 - - call diffusesmoothonestep(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & - tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & - ft(:,t),kt(:,:,t), im,p,m,r,p,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat(:,t),needeps) - - - + call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & + ft(:,t),kt(:,:,t), im,p,m,r,p,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat(:,t),needeps) end do end if diff --git a/src/smoothsimfast.f95 b/src/smoothsimfast.f95 index c87bc67..e89056e 100644 --- a/src/smoothsimfast.f95 +++ b/src/smoothsimfast.f95 @@ -6,7 +6,7 @@ subroutine smoothsimfast(yt, ymiss, timevar, zt, ht,tt, rtv,qt,a1, ft,kt,& logical, intent(in) :: needeps integer, intent(in) :: p, m, r,n,dt,jt - integer :: t, i,j + integer :: t, i integer, intent(in), dimension(n,p) :: ymiss integer, intent(in), dimension(5) :: timevar double precision, intent(in), dimension(n,p) :: yt @@ -24,33 +24,32 @@ subroutine smoothsimfast(yt, ymiss, timevar, zt, ht,tt, rtv,qt,a1, ft,kt,& double precision, dimension(m) :: at double precision, dimension(m,m) :: im double precision, intent(inout), dimension(m) :: rt0,rt1 - double precision, external :: ddot - - external dgemv, dger, dsymv - + double precision :: lik + lik = 0.0d0 at = a1 - !diffuse filtering begins - do t = 1, dt - 1 - call diffusefilteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& - tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),& - finf(:,t),kinf(:,:,t),p,m,p) - end do - - t = dt - call diffusefilteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& - tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),& - finf(:,t),kinf(:,:,t),p,m,jt) - !non-diffuse filtering begins - call filteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& - tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),p,m,jt) - - if(dt.LT.n) then - !Non-diffuse filtering continues from t=d+1, i=1 - do t = dt + 1, n - call filteronestepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& - tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),p,m,0) + if(dt .GT. 0) then + !diffuse filtering begins + do t = 1, dt - 1 + call dfilter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),& + finf(:,t),kinf(:,:,t),p,m,p,lik) end do + + t = dt + call dfilter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),& + finf(:,t),kinf(:,:,t),p,m,jt,lik) + !non-diffuse filtering begins + if(jt .LT. p) then + call filter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),p,m,jt,lik) + end if end if + !Non-diffuse filtering continues from t=d+1, i=1 + do t = dt + 1, n + call filter1stepnv(ymiss(t,:),yt(t,:),transpose(zt(:,:,(t-1)*timevar(1)+1)),& + tt(:,:,(t-1)*timevar(3)+1),at,vt(:,t),ft(:,t),kt(:,:,t),p,m,0,lik) + end do !smoothing begins @@ -63,24 +62,24 @@ subroutine smoothsimfast(yt, ymiss, timevar, zt, ht,tt, rtv,qt,a1, ft,kt,& rt0 = 0.0d0 do t = n, dt+1, -1 - call smoothonestep(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & ft(:,t),kt(:,:,t), im,p,m,r,1,rt0,etahat(:,t),epshat(:,t),needeps) end do - if(dt.GT.0) then - t=dt + if(dt .GT. 0) then + t = dt if(jt .LT. p) then - call smoothonestep(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + call smooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & ft(:,t),kt(:,:,t), im,p,m,r,jt+1,rt0,etahat(:,t),epshat(:,t),needeps) end if rt1 = 0.0d0 - call diffusesmoothonestep(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & ft(:,t),kt(:,:,t), im,p,m,r,jt,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat(:,t),needeps) - do t=(dt-1), 1, -1 - call diffusesmoothonestep(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & + do t = (dt - 1), 1, -1 + call dsmooth1step(ymiss(t,:), transpose(zt(:,:,(t-1)*timevar(1)+1)), ht(:,:,(t-1)*timevar(2)+1), & tt(:,:,(t-1)*timevar(3)+1), rtv(:,:,(t-1)*timevar(4)+1), qt(:,:,(t-1)*timevar(5)+1), vt(:,t), & ft(:,t),kt(:,:,t), im,p,m,r,p,rt0,rt1,finf(:,t),kinf(:,:,t),etahat(:,t),epshat(:,t),needeps) end do diff --git a/vignettes/KFAS-concordance.tex b/vignettes/KFAS-concordance.tex new file mode 100644 index 0000000..2c82c04 --- /dev/null +++ b/vignettes/KFAS-concordance.tex @@ -0,0 +1,3 @@ +\Sconcordance{concordance:KFAS.tex:KFAS.Rnw:% +1 34 1 49 0 1 6 387 1 4 0 8 1 4 0 41 1 4 0 16 1 17 0 59 1 30 0 % +11 1 21 0 7 1 9 0 8 1 9 0 13 1 4 0 13 1 9 0 263 1} diff --git a/vignettes/KFAS.log b/vignettes/KFAS.log new file mode 100644 index 0000000..e9ae050 --- /dev/null +++ b/vignettes/KFAS.log @@ -0,0 +1,759 @@ +This is pdfTeX, Version 3.1415926-2.5-1.40.14 (MiKTeX 2.9) (preloaded format=pdflatex 2014.10.27) 15 APR 2015 16:25 +entering extended mode +**KFAS.tex +(U:\MyPrograms\gitrepos\KFAS\vignettes\KFAS.tex +LaTeX2e <2011/06/27> +Babel and hyphenation patterns for english, afrikaans, ancientgreek, ar +abic, armenian, assamese, basque, bengali, bokmal, bulgarian, catalan, coptic, +croatian, czech, danish, dutch, esperanto, estonian, farsi, finnish, french, ga +lician, german, german-x-2013-05-26, greek, gujarati, hindi, hungarian, iceland +ic, indonesian, interlingua, irish, italian, kannada, kurmanji, latin, latvian, + lithuanian, malayalam, marathi, mongolian, mongolianlmc, monogreek, ngerman, n +german-x-2013-05-26, nynorsk, oriya, panjabi, pinyin, polish, portuguese, roman +ian, russian, sanskrit, serbian, slovak, slovenian, spanish, swedish, swissgerm +an, tamil, telugu, turkish, turkmen, ukenglish, ukrainian, uppersorbian, usengl +ishmax, welsh, loaded. +(C:/PROGRA~1/R/share/texmf/tex/latex\jss.cls +Document Class: jss 2013/09/05 2.3 jss class by Achim Zeileis +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\article.cls" +Document Class: article 2007/10/19 v1.4h Standard LaTeX document class +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\size11.clo" +File: size11.clo 2007/10/19 v1.4h Standard LaTeX file (size option) +) +\c@part=\count79 +\c@section=\count80 +\c@subsection=\count81 +\c@subsubsection=\count82 +\c@paragraph=\count83 +\c@subparagraph=\count84 +\c@figure=\count85 +\c@table=\count86 +\abovecaptionskip=\skip41 +\belowcaptionskip=\skip42 +\bibindent=\dimen102 +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\graphics\graphicx.sty" +Package: graphicx 1999/02/16 v1.0f Enhanced LaTeX Graphics (DPC,SPQR) + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\graphics\keyval.sty" +Package: keyval 1999/03/16 v1.13 key=value parser (DPC) +\KV@toks@=\toks14 +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\graphics\graphics.sty" +Package: graphics 2009/02/05 v1.0o Standard LaTeX Graphics (DPC,SPQR) + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\graphics\trig.sty" +Package: trig 1999/03/16 v1.09 sin cos tan (DPC) +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\00miktex\graphics.cfg" +File: graphics.cfg 2007/01/18 v1.5 graphics configuration of teTeX/TeXLive +) +Package graphics Info: Driver file: pdftex.def on input line 91. + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\pdftex-def\pdftex.def" +File: pdftex.def 2011/05/27 v0.06d Graphics/color for pdfTeX + +("C:\Program Files (x86)\MiKTeX 2.9\tex\generic\oberdiek\infwarerr.sty" +Package: infwarerr 2010/04/08 v1.3 Providing info/warning/error messages (HO) +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\generic\oberdiek\ltxcmds.sty" +Package: ltxcmds 2011/11/09 v1.22 LaTeX kernel commands for general use (HO) +) +\Gread@gobject=\count87 +)) +\Gin@req@height=\dimen103 +\Gin@req@width=\dimen104 +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\graphics\color.sty" +Package: color 2005/11/14 v1.0j Standard LaTeX Color (DPC) + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\00miktex\color.cfg" +File: color.cfg 2007/01/18 v1.5 color configuration of teTeX/TeXLive +) +Package color Info: Driver file: pdftex.def on input line 130. +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\ae\ae.sty" +Package: ae 2001/02/12 1.3 Almost European Computer Modern + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\fontenc.sty" +Package: fontenc 2005/09/27 v1.99g Standard LaTeX package + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\t1enc.def" +File: t1enc.def 2005/09/27 v1.99g Standard LaTeX file +LaTeX Font Info: Redeclaring font encoding T1 on input line 43. +) +LaTeX Font Info: Try loading font information for T1+aer on input line 100. + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\ae\t1aer.fd" +File: t1aer.fd 1997/11/16 Font definitions for T1/aer. +))) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\fancyvrb\fancyvrb.sty" +Package: fancyvrb 2008/02/07 + +Style option: `fancyvrb' v2.7a, with DG/SPQR fixes, and firstline=lastline fix +<2008/02/07> (tvz) +\FV@CodeLineNo=\count88 +\FV@InFile=\read1 +\FV@TabBox=\box26 +\c@FancyVerbLine=\count89 +\FV@StepNumber=\count90 +\FV@OutFile=\write3 +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\fontenc.sty" +Package: fontenc 2005/09/27 v1.99g Standard LaTeX package + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\t1enc.def" +File: t1enc.def 2005/09/27 v1.99g Standard LaTeX file +LaTeX Font Info: Redeclaring font encoding T1 on input line 43. +)) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\upquote\upquote.sty" +Package: upquote 2012/04/19 v1.3 upright-quote and grave-accent glyphs in verba +tim + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\textcomp.sty" +Package: textcomp 2005/09/27 v1.99g Standard LaTeX package +Package textcomp Info: Sub-encoding information: +(textcomp) 5 = only ISO-Adobe without \textcurrency +(textcomp) 4 = 5 + \texteuro +(textcomp) 3 = 4 + \textohm +(textcomp) 2 = 3 + \textestimated + \textcurrency +(textcomp) 1 = TS1 - \textcircled - \t +(textcomp) 0 = TS1 (full) +(textcomp) Font families with sub-encoding setting implement +(textcomp) only a restricted character set as indicated. +(textcomp) Family '?' is the default used for unknown fonts. +(textcomp) See the documentation for details. +Package textcomp Info: Setting ? sub-encoding to TS1/1 on input line 71. + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\ts1enc.def" +File: ts1enc.def 2001/06/05 v3.0e (jk/car/fm) Standard LaTeX file +) +LaTeX Info: Redefining \oldstylenums on input line 266. +Package textcomp Info: Setting cmr sub-encoding to TS1/0 on input line 281. +Package textcomp Info: Setting cmss sub-encoding to TS1/0 on input line 282. +Package textcomp Info: Setting cmtt sub-encoding to TS1/0 on input line 283. +Package textcomp Info: Setting cmvtt sub-encoding to TS1/0 on input line 284. +Package textcomp Info: Setting cmbr sub-encoding to TS1/0 on input line 285. +Package textcomp Info: Setting cmtl sub-encoding to TS1/0 on input line 286. +Package textcomp Info: Setting ccr sub-encoding to TS1/0 on input line 287. +Package textcomp Info: Setting ptm sub-encoding to TS1/4 on input line 288. +Package textcomp Info: Setting pcr sub-encoding to TS1/4 on input line 289. +Package textcomp Info: Setting phv sub-encoding to TS1/4 on input line 290. +Package textcomp Info: Setting ppl sub-encoding to TS1/3 on input line 291. +Package textcomp Info: Setting pag sub-encoding to TS1/4 on input line 292. +Package textcomp Info: Setting pbk sub-encoding to TS1/4 on input line 293. +Package textcomp Info: Setting pnc sub-encoding to TS1/4 on input line 294. +Package textcomp Info: Setting pzc sub-encoding to TS1/4 on input line 295. +Package textcomp Info: Setting bch sub-encoding to TS1/4 on input line 296. +Package textcomp Info: Setting put sub-encoding to TS1/5 on input line 297. +Package textcomp Info: Setting uag sub-encoding to TS1/5 on input line 298. +Package textcomp Info: Setting ugq sub-encoding to TS1/5 on input line 299. +Package textcomp Info: Setting ul8 sub-encoding to TS1/4 on input line 300. +Package textcomp Info: Setting ul9 sub-encoding to TS1/4 on input line 301. +Package textcomp Info: Setting augie sub-encoding to TS1/5 on input line 302. +Package textcomp Info: Setting dayrom sub-encoding to TS1/3 on input line 303. +Package textcomp Info: Setting dayroms sub-encoding to TS1/3 on input line 304. + +Package textcomp Info: Setting pxr sub-encoding to TS1/0 on input line 305. +Package textcomp Info: Setting pxss sub-encoding to TS1/0 on input line 306. +Package textcomp Info: Setting pxtt sub-encoding to TS1/0 on input line 307. +Package textcomp Info: Setting txr sub-encoding to TS1/0 on input line 308. +Package textcomp Info: Setting txss sub-encoding to TS1/0 on input line 309. +Package textcomp Info: Setting txtt sub-encoding to TS1/0 on input line 310. +Package textcomp Info: Setting lmr sub-encoding to TS1/0 on input line 311. +Package textcomp Info: Setting lmdh sub-encoding to TS1/0 on input line 312. +Package textcomp Info: Setting lmss sub-encoding to TS1/0 on input line 313. +Package textcomp Info: Setting lmssq sub-encoding to TS1/0 on input line 314. +Package textcomp Info: Setting lmvtt sub-encoding to TS1/0 on input line 315. +Package textcomp Info: Setting qhv sub-encoding to TS1/0 on input line 316. +Package textcomp Info: Setting qag sub-encoding to TS1/0 on input line 317. +Package textcomp Info: Setting qbk sub-encoding to TS1/0 on input line 318. +Package textcomp Info: Setting qcr sub-encoding to TS1/0 on input line 319. +Package textcomp Info: Setting qcs sub-encoding to TS1/0 on input line 320. +Package textcomp Info: Setting qpl sub-encoding to TS1/0 on input line 321. +Package textcomp Info: Setting qtm sub-encoding to TS1/0 on input line 322. +Package textcomp Info: Setting qzc sub-encoding to TS1/0 on input line 323. +Package textcomp Info: Setting qhvc sub-encoding to TS1/0 on input line 324. +Package textcomp Info: Setting futs sub-encoding to TS1/4 on input line 325. +Package textcomp Info: Setting futx sub-encoding to TS1/4 on input line 326. +Package textcomp Info: Setting futj sub-encoding to TS1/4 on input line 327. +Package textcomp Info: Setting hlh sub-encoding to TS1/3 on input line 328. +Package textcomp Info: Setting hls sub-encoding to TS1/3 on input line 329. +Package textcomp Info: Setting hlst sub-encoding to TS1/3 on input line 330. +Package textcomp Info: Setting hlct sub-encoding to TS1/5 on input line 331. +Package textcomp Info: Setting hlx sub-encoding to TS1/5 on input line 332. +Package textcomp Info: Setting hlce sub-encoding to TS1/5 on input line 333. +Package textcomp Info: Setting hlcn sub-encoding to TS1/5 on input line 334. +Package textcomp Info: Setting hlcw sub-encoding to TS1/5 on input line 335. +Package textcomp Info: Setting hlcf sub-encoding to TS1/5 on input line 336. +Package textcomp Info: Setting pplx sub-encoding to TS1/3 on input line 337. +Package textcomp Info: Setting pplj sub-encoding to TS1/3 on input line 338. +Package textcomp Info: Setting ptmx sub-encoding to TS1/4 on input line 339. +Package textcomp Info: Setting ptmj sub-encoding to TS1/4 on input line 340. +)) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\natbib\natbib.sty" +Package: natbib 2010/09/13 8.31b (PWD, AO) +\bibhang=\skip43 +\bibsep=\skip44 +LaTeX Info: Redefining \cite on input line 694. +\c@NAT@ctr=\count91 +) +\footerskip=\skip45 + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\hyperref\hyperref.sty" +Package: hyperref 2012/11/06 v6.83m Hypertext links for LaTeX + +("C:\Program Files (x86)\MiKTeX 2.9\tex\generic\oberdiek\hobsub-hyperref.sty" +Package: hobsub-hyperref 2012/04/25 v1.12 Bundle oberdiek, subset hyperref (HO) + + +("C:\Program Files (x86)\MiKTeX 2.9\tex\generic\oberdiek\hobsub-generic.sty" +Package: hobsub-generic 2012/04/25 v1.12 Bundle oberdiek, subset generic (HO) +Package: hobsub 2012/04/25 v1.12 Construct package bundles (HO) +Package hobsub Info: Skipping package `infwarerr' (already loaded). +Package hobsub Info: Skipping package `ltxcmds' (already loaded). +Package: ifluatex 2010/03/01 v1.3 Provides the ifluatex switch (HO) +Package ifluatex Info: LuaTeX not detected. +Package: ifvtex 2010/03/01 v1.5 Detect VTeX and its facilities (HO) +Package ifvtex Info: VTeX not detected. +Package: intcalc 2007/09/27 v1.1 Expandable calculations with integers (HO) +Package: ifpdf 2011/01/30 v2.3 Provides the ifpdf switch (HO) +Package ifpdf Info: pdfTeX in PDF mode is detected. +Package: etexcmds 2011/02/16 v1.5 Avoid name clashes with e-TeX commands (HO) +Package etexcmds Info: Could not find \expanded. +(etexcmds) That can mean that you are not using pdfTeX 1.50 or +(etexcmds) that some package has redefined \expanded. +(etexcmds) In the latter case, load this package earlier. +Package: kvsetkeys 2012/04/25 v1.16 Key value parser (HO) +Package: kvdefinekeys 2011/04/07 v1.3 Define keys (HO) +Package: pdftexcmds 2011/11/29 v0.20 Utility functions of pdfTeX for LuaTeX (HO +) +Package pdftexcmds Info: LuaTeX not detected. +Package pdftexcmds Info: \pdf@primitive is available. +Package pdftexcmds Info: \pdf@ifprimitive is available. +Package pdftexcmds Info: \pdfdraftmode found. +Package: pdfescape 2011/11/25 v1.13 Implements pdfTeX's escape features (HO) +Package: bigintcalc 2012/04/08 v1.3 Expandable calculations on big integers (HO +) +Package: bitset 2011/01/30 v1.1 Handle bit-vector datatype (HO) +Package: uniquecounter 2011/01/30 v1.2 Provide unlimited unique counter (HO) +) +Package hobsub Info: Skipping package `hobsub' (already loaded). +Package: letltxmacro 2010/09/02 v1.4 Let assignment for LaTeX macros (HO) +Package: hopatch 2011/06/24 v1.1 Wrapper for package hooks (HO) +Package: xcolor-patch 2011/01/30 xcolor patch +Package: atveryend 2011/06/30 v1.8 Hooks at the very end of document (HO) +Package atveryend Info: \enddocument detected (standard20110627). +Package: atbegshi 2011/10/05 v1.16 At begin shipout hook (HO) +Package: refcount 2011/10/16 v3.4 Data extraction from label references (HO) +Package: hycolor 2011/01/30 v1.7 Color options for hyperref/bookmark (HO) +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\generic\ifxetex\ifxetex.sty" +Package: ifxetex 2010/09/12 v0.6 Provides ifxetex conditional +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\oberdiek\auxhook.sty" +Package: auxhook 2011/03/04 v1.3 Hooks for auxiliary files (HO) +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\oberdiek\kvoptions.sty" +Package: kvoptions 2011/06/30 v3.11 Key value format for package options (HO) +) +\@linkdim=\dimen105 +\Hy@linkcounter=\count92 +\Hy@pagecounter=\count93 + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\hyperref\pd1enc.def" +File: pd1enc.def 2012/11/06 v6.83m Hyperref: PDFDocEncoding definition (HO) +) +\Hy@SavedSpaceFactor=\count94 + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\00miktex\hyperref.cfg" +File: hyperref.cfg 2002/06/06 v1.2 hyperref configuration of TeXLive +) +Package hyperref Info: Hyper figures OFF on input line 4443. +Package hyperref Info: Link nesting OFF on input line 4448. +Package hyperref Info: Hyper index ON on input line 4451. +Package hyperref Info: Plain pages OFF on input line 4458. +Package hyperref Info: Backreferencing OFF on input line 4463. +Package hyperref Info: Implicit mode ON; LaTeX internals redefined. +Package hyperref Info: Bookmarks ON on input line 4688. +\c@Hy@tempcnt=\count95 + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\url\url.sty" +\Urlmuskip=\muskip10 +Package: url 2013/09/16 ver 3.4 Verb mode for urls, etc. +) +LaTeX Info: Redefining \url on input line 5041. +\XeTeXLinkMargin=\dimen106 +\Fld@menulength=\count96 +\Field@Width=\dimen107 +\Fld@charsize=\dimen108 +Package hyperref Info: Hyper figures OFF on input line 6295. +Package hyperref Info: Link nesting OFF on input line 6300. +Package hyperref Info: Hyper index ON on input line 6303. +Package hyperref Info: backreferencing OFF on input line 6310. +Package hyperref Info: Link coloring OFF on input line 6315. +Package hyperref Info: Link coloring with OCG OFF on input line 6320. +Package hyperref Info: PDF/A mode OFF on input line 6325. +LaTeX Info: Redefining \ref on input line 6365. +LaTeX Info: Redefining \pageref on input line 6369. +\Hy@abspage=\count97 +\c@Item=\count98 +\c@Hfootnote=\count99 +) + +Package hyperref Message: Driver (autodetected): hpdftex. + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\hyperref\hpdftex.def" +File: hpdftex.def 2012/11/06 v6.83m Hyperref driver for pdfTeX +\Fld@listcount=\count100 +\c@bookmark@seq@number=\count101 + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\oberdiek\rerunfilecheck.sty" +Package: rerunfilecheck 2011/04/15 v1.7 Rerun checks for auxiliary files (HO) +Package uniquecounter Info: New unique counter `rerunfilecheck' on input line 2 +82. +) +\Hy@SectionHShift=\skip46 +) +\preXLskip=\skip47 +\preLskip=\skip48 +\preMskip=\skip49 +\preSskip=\skip50 +\postMskip=\skip51 +\postSskip=\skip52 + + +Package hyperref Warning: Option `hyperindex' has already been used, +(hyperref) setting the option has no effect on input line 444. + +Package hyperref Info: Option `colorlinks' set `true' on input line 444. +Package hyperref Info: Option `linktocpage' set `true' on input line 444. +Package hyperref Info: Option `plainpages' set `false' on input line 444. +) ("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\framed\framed.sty" +Package: framed 2011/10/22 v 0.96: framed or shaded text with page breaks +\OuterFrameSep=\skip53 +\fb@frw=\dimen109 +\fb@frh=\dimen110 +\FrameRule=\dimen111 +\FrameSep=\dimen112 +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\alltt.sty" +Package: alltt 1997/06/16 v2.0g defines alltt environment +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\amsmath\amsmath.sty" +Package: amsmath 2013/01/14 v2.14 AMS math features +\@mathmargin=\skip54 + +For additional information on amsmath, use the `?' option. +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\amsmath\amstext.sty" +Package: amstext 2000/06/29 v2.01 + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\amsmath\amsgen.sty" +File: amsgen.sty 1999/11/30 v2.0 +\@emptytoks=\toks15 +\ex@=\dimen113 +)) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\amsmath\amsbsy.sty" +Package: amsbsy 1999/11/29 v1.2d +\pmbraise@=\dimen114 +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\amsmath\amsopn.sty" +Package: amsopn 1999/12/14 v2.01 operator names +) +\inf@bad=\count102 +LaTeX Info: Redefining \frac on input line 210. +\uproot@=\count103 +\leftroot@=\count104 +LaTeX Info: Redefining \overline on input line 306. +\classnum@=\count105 +\DOTSCASE@=\count106 +LaTeX Info: Redefining \ldots on input line 378. +LaTeX Info: Redefining \dots on input line 381. +LaTeX Info: Redefining \cdots on input line 466. +\Mathstrutbox@=\box27 +\strutbox@=\box28 +\big@size=\dimen115 +LaTeX Font Info: Redeclaring font encoding OML on input line 566. +LaTeX Font Info: Redeclaring font encoding OMS on input line 567. +\macc@depth=\count107 +\c@MaxMatrixCols=\count108 +\dotsspace@=\muskip11 +\c@parentequation=\count109 +\dspbrk@lvl=\count110 +\tag@help=\toks16 +\row@=\count111 +\column@=\count112 +\maxfields@=\count113 +\andhelp@=\toks17 +\eqnshift@=\dimen116 +\alignsep@=\dimen117 +\tagshift@=\dimen118 +\tagwidth@=\dimen119 +\totwidth@=\dimen120 +\lineht@=\dimen121 +\@envbody=\toks18 +\multlinegap=\skip55 +\multlinetaggap=\skip56 +\mathdisplay@stack=\toks19 +LaTeX Info: Redefining \[ on input line 2665. +LaTeX Info: Redefining \] on input line 2666. +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\ulem\ulem.sty" +\UL@box=\box29 +\UL@hyphenbox=\box30 +\UL@skip=\skip57 +\UL@hook=\toks20 +\UL@height=\dimen122 +\UL@pe=\count114 +\UL@pixel=\dimen123 +\ULC@box=\box31 +Package: ulem 2012/05/18 +\ULdepth=\dimen124 +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\inputenc.sty" +Package: inputenc 2008/03/30 v1.1d Input encoding file +\inpenc@prehook=\toks21 +\inpenc@posthook=\toks22 + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\latin1.def" +File: latin1.def 2008/03/30 v1.1d Input encoding file +)) +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\fontenc.sty" +Package: fontenc 2005/09/27 v1.99g Standard LaTeX package + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\t1enc.def" +File: t1enc.def 2005/09/27 v1.99g Standard LaTeX file +LaTeX Font Info: Redeclaring font encoding T1 on input line 43. +)) +("C:\Program Files (x86)\MiKTeX 2.9\tex\generic\babel\babel.sty" +Package: babel 2008/07/08 v3.8m The Babel package + +************************************* +* Local config file bblopts.cfg used +* +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\00miktex\bblopts.cfg" +File: bblopts.cfg 2006/07/31 v1.0 MiKTeX 'babel' configuration +) +("C:\Program Files (x86)\MiKTeX 2.9\tex\generic\babel\english.ldf" +Language: english 2005/03/30 v3.3o English support from the babel system + +("C:\Program Files (x86)\MiKTeX 2.9\tex\generic\babel\babel.def" +File: babel.def 2008/07/08 v3.8m Babel common definitions +\babel@savecnt=\count115 +\U@D=\dimen125 +) +\l@canadian = a dialect from \language\l@american +\l@australian = a dialect from \language\l@british +\l@newzealand = a dialect from \language\l@british +)) +(U:\MyPrograms\gitrepos\KFAS\vignettes\KFAS.aux) +LaTeX Font Info: Checking defaults for OML/cmm/m/it on input line 82. +LaTeX Font Info: ... okay on input line 82. +LaTeX Font Info: Checking defaults for T1/cmr/m/n on input line 82. +LaTeX Font Info: ... okay on input line 82. +LaTeX Font Info: Checking defaults for OT1/cmr/m/n on input line 82. +LaTeX Font Info: ... okay on input line 82. +LaTeX Font Info: Checking defaults for OMS/cmsy/m/n on input line 82. +LaTeX Font Info: ... okay on input line 82. +LaTeX Font Info: Checking defaults for OMX/cmex/m/n on input line 82. +LaTeX Font Info: ... okay on input line 82. +LaTeX Font Info: Checking defaults for U/cmr/m/n on input line 82. +LaTeX Font Info: ... okay on input line 82. +LaTeX Font Info: Checking defaults for TS1/cmr/m/n on input line 82. +LaTeX Font Info: Try loading font information for TS1+cmr on input line 82. + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\base\ts1cmr.fd" +File: ts1cmr.fd 1999/05/25 v2.5h Standard LaTeX font definitions +) +LaTeX Font Info: ... okay on input line 82. +LaTeX Font Info: Checking defaults for PD1/pdf/m/n on input line 82. +LaTeX Font Info: ... okay on input line 82. + +("C:\Program Files (x86)\MiKTeX 2.9\tex\context\base\supp-pdf.mkii" +[Loading MPS to PDF converter (version 2006.09.02).] +\scratchcounter=\count116 +\scratchdimen=\dimen126 +\scratchbox=\box32 +\nofMPsegments=\count117 +\nofMParguments=\count118 +\everyMPshowfont=\toks23 +\MPscratchCnt=\count119 +\MPscratchDim=\dimen127 +\MPnumerator=\count120 +\makeMPintoPDFobject=\count121 +\everyMPtoPDFconversion=\toks24 +) +\AtBeginShipoutBox=\box33 +Package hyperref Info: Link coloring ON on input line 82. + ("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\hyperref\nameref.sty" +Package: nameref 2012/10/27 v2.43 Cross-referencing by name of section + +("C:\Program Files (x86)\MiKTeX 2.9\tex\generic\oberdiek\gettitlestring.sty" +Package: gettitlestring 2010/12/03 v1.4 Cleanup title references (HO) +) +\c@section@level=\count122 +) +LaTeX Info: Redefining \ref on input line 82. +LaTeX Info: Redefining \pageref on input line 82. +LaTeX Info: Redefining \nameref on input line 82. + +(U:\MyPrograms\gitrepos\KFAS\vignettes\KFAS.out) +(U:\MyPrograms\gitrepos\KFAS\vignettes\KFAS.out) +\@outlinefile=\write4 +LaTeX Font Info: Try loading font information for T1+aess on input line 82. + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\ae\t1aess.fd" +File: t1aess.fd 1997/11/16 Font definitions for T1/aess. +) [1 + +{C:/localtexmf/pdftex/config/pdftex.map}] +LaTeX Font Info: Try loading font information for T1+aett on input line 115. + + +("C:\Program Files (x86)\MiKTeX 2.9\tex\latex\ae\t1aett.fd" +File: t1aett.fd 1997/11/16 Font definitions for T1/aett. +) +Overfull \hbox (5.47499pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 2 [] + [] + +[2] +Overfull \hbox (5.47499pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 3 + [] + +[3] +Overfull \hbox (5.47499pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 4 [] + [] + +[4] +Overfull \hbox (5.47499pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 5 + [] + +[5] +Overfull \hbox (5.47499pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 6 [] + [] + +[6] +Overfull \hbox (5.47499pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 7 + [] + +[7] +Overfull \hbox (5.47499pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 8 [] + [] + +[8] +Underfull \vbox (badness 10000) has occurred while \output is active [] + + +Overfull \hbox (5.47499pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 9 + [] + +[9] +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 10 [] + [] + +[10] +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 11 + [] + +[11] +Package color Info: Redefining color shadecolor on input line 465. + + +LaTeX Font Warning: Font shape `T1/aett/bx/n' undefined +(Font) using `T1/aett/m/n' instead on input line 467. + +Package color Info: Redefining color shadecolor on input line 480. + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 12 [] + [] + +[12] +Package color Info: Redefining color shadecolor on input line 502. + +Overfull \hbox (1.04945pt too wide) in paragraph at lines 518--518 +[][]\T1/aett/m/it/10.95 # Com-mon fixed ef-fect part and dis-tinct ran-dom ef-f +ect parts for each "se-ries"[][] + [] + + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 13 + [] + +[13] +Package color Info: Redefining color shadecolor on input line 541. +
+File: figure/alcoholPlot1.pdf Graphic file (type pdf) + + +Package pdftex.def Info: figure/alcoholPlot1.pdf used on input line 561. +(pdftex.def) Requested size: 442.65375pt x 442.66037pt. + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 14 [] + [] + +[14] +Package color Info: Redefining color shadecolor on input line 586. + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 15 + [] + +[15 ] +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 16 [] + [] + +[16] +Package color Info: Redefining color shadecolor on input line 663. +Package color Info: Redefining color shadecolor on input line 691. + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 17 + [] + +[17]
+File: figure/states.pdf Graphic file (type pdf) + +Package pdftex.def Info: figure/states.pdf used on input line 698. +(pdftex.def) Requested size: 442.65375pt x 442.66037pt. +Package color Info: Redefining color shadecolor on input line 707. + +Underfull \vbox (badness 1308) has occurred while \output is active [] + + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 18 [] + [] + +[18 ] +
+File: figure/diagnostics1.pdf Graphic file (type pdf) + + +Package pdftex.def Info: figure/diagnostics1.pdf used on input line 715. +(pdftex.def) Requested size: 442.65375pt x 442.66037pt. + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 19 + [] + +[19 ] +Package color Info: Redefining color shadecolor on input line 726. +Package color Info: Redefining color shadecolor on input line 738. + +
+File: figure/predictplot.pdf Graphic file (type pdf) + + +Package pdftex.def Info: figure/predictplot.pdf used on input line 754. +(pdftex.def) Requested size: 442.65375pt x 442.66037pt. + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 20 [] + [] + +[20] +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 21 + [] + +[21 ] +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 22 [] + [] + +[22 + +] +Overfull \hbox (26.85487pt too wide) detected at line 964 +[] [] + [] + + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 23 + [] + +[23] (U:\MyPrograms\gitrepos\KFAS\vignettes\KFAS.bbl +LaTeX Font Info: Font shape `T1/aer/b/it' in size <10.95> not available +(Font) Font shape `T1/aer/bx/it' tried instead on input line 19. + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 24 [] + [] + +[24] +LaTeX Font Info: Font shape `T1/aess/m/it' in size <10.95> not available +(Font) Font shape `T1/aess/m/sl' tried instead on input line 33. + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +[] \T1/aer/m/n/10.95 25 + [] + +[25]) +Package atveryend Info: Empty hook `BeforeClearDocument' on input line 1020. + +Overfull \hbox (10.94998pt too wide) has occurred while \output is active +\T1/aer/m/n/10.95 26 [] + [] + +[26] +Package atveryend Info: Empty hook `AfterLastShipout' on input line 1020. + (U:\MyPrograms\gitrepos\KFAS\vignettes\KFAS.aux) +Package atveryend Info: Executing hook `AtVeryEndDocument' on input line 1020. +Package atveryend Info: Executing hook `AtEndAfterFileList' on input line 1020. + +Package rerunfilecheck Info: File `KFAS.out' has not changed. +(rerunfilecheck) Checksum: 9AE816EC4F57E52B7D1D45A84CF8983C;2035. + + +LaTeX Font Warning: Some font shapes were not available, defaults substituted. + +Package atveryend Info: Empty hook `AtVeryVeryEnd' on input line 1020. + ) +Here is how much of TeX's memory you used: + 7969 strings out of 493921 + 114199 string characters out of 3144866 + 203099 words of memory out of 3000000 + 10994 multiletter control sequences out of 15000+200000 + 52998 words of font info for 95 fonts, out of 3000000 for 9000 + 841 hyphenation exceptions out of 8191 + 36i,10n,34p,1482b,466s stack positions out of 5000i,500n,10000p,200000b,50000s +< +C:/Program Files (x86)/MiKTeX 2.9/fonts/type1/public/amsfonts/cm/cmbx10.pfb> +Output written on KFAS.pdf (26 pages, 418572 bytes). +PDF statistics: + 483 PDF objects out of 1000 (max. 8388607) + 99 named destinations out of 1000 (max. 500000) + 157 words of extra memory for PDF output out of 10000 (max. 10000000) + diff --git a/vignettes/KFAS.pdf b/vignettes/KFAS.pdf new file mode 100644 index 0000000..c58ab13 Binary files /dev/null and b/vignettes/KFAS.pdf differ diff --git a/vignettes/KFAS.synctex.gz b/vignettes/KFAS.synctex.gz new file mode 100644 index 0000000..2bec376 Binary files /dev/null and b/vignettes/KFAS.synctex.gz differ diff --git a/vignettes/KFAS.tex b/vignettes/KFAS.tex new file mode 100644 index 0000000..5b8f490 --- /dev/null +++ b/vignettes/KFAS.tex @@ -0,0 +1,1021 @@ +%\VignetteEngine{knitr::knitr} +%\VignetteDepends{lme4} +%\VignetteIndexEntry{KFAS: Exponential Family State Space Models in R} +\documentclass[nojss,article]{jss}\usepackage[]{graphicx}\usepackage[]{color} +%% maxwidth is the original width if it is less than linewidth +%% otherwise use linewidth (to make sure the graphics do not exceed the margin) +\makeatletter +\def\maxwidth{ % + \ifdim\Gin@nat@width>\linewidth + \linewidth + \else + \Gin@nat@width + \fi +} +\makeatother + +\definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345} +\newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}% +\newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}% +\newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}% +\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}% +\newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}% +\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}% +\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}% +\newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}% +\newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}% + +\usepackage{framed} +\makeatletter +\newenvironment{kframe}{% + \def\at@end@of@kframe{}% + \ifinner\ifhmode% + \def\at@end@of@kframe{\end{minipage}}% + \begin{minipage}{\columnwidth}% + \fi\fi% + \def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep + \colorbox{shadecolor}{##1}\hskip-\fboxsep + % There is no \\@totalrightmargin, so: + \hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}% + \MakeFramed {\advance\hsize-\width + \@totalleftmargin\z@ \linewidth\hsize + \@setminipage}}% + {\par\unskip\endMakeFramed% + \at@end@of@kframe} +\makeatother + +\definecolor{shadecolor}{rgb}{.97, .97, .97} +\definecolor{messagecolor}{rgb}{0, 0, 0} +\definecolor{warningcolor}{rgb}{1, 0, 1} +\definecolor{errorcolor}{rgb}{1, 0, 0} +\newenvironment{knitrout}{}{} % an empty environment to be redefined in TeX + +\usepackage{alltt} +\usepackage{amsmath} +\usepackage{natbib} +\usepackage[normalem]{ulem} +\usepackage[latin1]{inputenc} +\usepackage[T1]{fontenc} +\usepackage[english]{babel} +\author{Jouni Helske\\ University of Jyv\"askyl\"a} +\title{\pkg{KFAS}: Exponential family state space models in R} + +\Plainauthor{Jouni Helske} %% comma-separated +\Plaintitle{KFAS: Exponential Family State Space Models in R} %% without formatting +\Shorttitle{KFAS: Exponential Family State Space Models in R} %% a short title (if necessary) + +\Abstract{ +State space modelling is an efficient and flexible method for statistical inference of broad class of time series and other data. This paper describes an R package \pkg{KFAS} for modelling state space models with the observations from exponential family, namely Gaussian, Poisson, binomial, negative binomial and gamma distributions. After introducing the basic theory behind the state space modelling, an illustrative example is provided, and finally a short comparison to alternative modelling framework is presented. +} + +\Keywords{state space models, time series, dynamic linear models, \pkg{KFAS}, \proglang{R}, exponential family} +\Plainkeywords{state space models, exponential family, dynamic linear models, time series, KFAS, R} %% without formatting + +\Address{ +Jouni Helske\\ +University of Jyv\"askyl\"a\\ +Department of Mathematics and Statistics\\ +40014 Jyv\"askyl\"a, Finland\\ +E-mail: \email{Jouni.Helske@jyu.fi} +} +\IfFileExists{upquote.sty}{\usepackage{upquote}}{} +\begin{document} + + + +Modification of manuscript submitted to \emph{Journal of Statistical Software} + + +\section{Introduction} + +State space models offer an unified framework for modelling several types of time series and other data. Structural time series, ARIMA models, simple regression, generalized linear mixed models, and cubic spline smoothing are just some examples of the statistical models which can be represented as a state space model. One of the simplest classes of state space models are linear Gaussian state space models (also known as dynamic linear models), which are analytically tractable, and are therefore often used in many fields of science. + +\citet{Petris2010} and \citet{Tusell2010} introduce and review some of the contributed packages available at Comprehensive R Archive Network (CRAN) for \proglang{R} \citep{R} for Gaussian state space modelling. Since then, several new additions have emerged in CRAN. Most of these packages use one or multiple packages reviewed in \citet{Tusell2010} for filtering and smoothing, and add new user interface and functionality for certain type of models. For example, package \pkg{rucm} \citep{rucm} is focused on structural time series, while \pkg{dlmodeler} \citep{dlmodeler} provides unified interface compatible with multiple packages, and \pkg{MARSS} \citep{MARSS,MARSSarticle} provides functions for a maximum likelihood estimation of large class of Gaussian state space models via the EM-algorithm. + +One of the packages reviewed in the aforementioned papers is \pkg{KFAS} (Kalman Filtering And Smoothing) which, in addition, of modelling the general linear Gaussian state space models, it can also be used in cases where the observations are from other exponential family models, namely binomial, Poisson, negative binomial, and Gamma models. To my knowledge, \pkg{KFAS} is currently the only package which supports non-Gaussian models in a classical state space modelling framework. Package \pkg{sspir} \citep{sspir} included similar capabilities, but it was removed from CRAN in 2013. In addition to \pkg{KFAS}, there are at least two packages which are designed for more broader class of problems, but can also be used for exponential family state space modelling. Package \pkg{pomp} \citep{pomp} offers functions for inference of state space models with non-Gaussian and non-linear observation and state equations by particle filtering and related methods. Another package suitable for state space modelling is \pkg{INLA} \citep{inla} (not available on CRAN), which can be used for Bayesian analysis via integrated nested Laplace approximation technique. Although it is often used in a spatial modelling via Gaussian random fields, it can also be used for certain temporal state space models where the state transitions are Gaussian. + +After the papers by \citet{Petris2010} and \citet{Tusell2010}, \pkg{KFAS} has been completely rewritten. Package is now much more user friendly due to the use of \proglang{R}'s symbolic formulas in model definition. The non-Gaussian modelling, which was somewhat experimental in old versions of \pkg{KFAS}, is now fully functional supporting multivariate models with different distributions. Many other features have also been added (such as methods for computing model residuals), performance of the main functions have improved and in the process several bugs have been also fixed. + +In this paper I first introduce the basic theory related to state space modelling, and then proceed to show main aspects of \pkg{KFAS} in more detail, illustrate its functionality by applying it to real life dataset \sout{and finally make short comparison between \pkg{INLA} and \pkg{KFAS}} (comparison to \pkg{INLA} is removed from vignette as \pkg{INLA} is not available at CRAN). + +\section{Gaussian state space model} + +The theory behind \pkg{KFAS} is mostly based on \citet{DK2012} and related articles by the same authors, and therefore the basic notation is nearly identical with the one used by Durbin and Koopman. +For linear Gaussian state space model with continuous states and discrete time intervals $t=1,\ldots,n$ we have +\begin{equation}\label{ssgeneral} +\begin{aligned} +y_t &= Z_t\alpha_t + \epsilon_t, \quad \textrm{(observation equation)}\\ +\alpha_{t+1} &= T_t\alpha_t + R_t\eta_{t}, \quad \textrm{(state equation)} +\end{aligned} +\end{equation} +where $\epsilon_{t} \sim N(0,H_t), \eta_{t} \sim N(0,Q_t)$ and $\alpha_1 \sim +N(a_1,P_1)$ independently of each other. We assume that $y_t$ is $p\times 1$, +$\alpha_{t+1}$ is $m \times 1$ and $\eta_{t}$ is $k \times 1$ vector. We also denote $\alpha = (\alpha^\top_1,\ldots,\alpha^\top_n)^\top$ and similarly $y = (y^\top_1,\ldots,y^\top_n)^\top$. + +Here $y_t$ contains the observations at time $t$, whereas $\alpha_t$ is a vector of latent state process at time point $t$. The system matrices $Z_t$, $T_t$, and $R_t$, together with the covariance matrices $H_t$ and $Q_t$ depend on the particular model definition, and are often time invariant, i.e., do not depend on $t$. Usually at least some of these matrices contain unknown parameters which need to be estimated. In \pkg{KFAS} one defines the model with the function \code{SSModel}. Function \code{SSModel} only builds the model, and does not perform estimation of unknown parameters, which differs from functions like \code{lm} which builds and estimates the model with one command. + +The main goal of the state space modelling is to gain knowledge of the latent states $\alpha$ given the observations $y$. This is achieved by using two important recursive algorithms, Kalman filtering and smoothing. From Kalman filtering algorithm we obtain the one step ahead predictions and the prediction errors + +\begin{equation*} +\begin{aligned} +a_{t+1} &= \E(\alpha_{t+1}|y_t,\ldots,y_1), \\ +v_{t} &= y_t - Z_t a_t, +\end{aligned} +\end{equation*} +and the related covariance matrices +\begin{equation*} +\begin{aligned} +P_{t+1} &= \VAR(\alpha_{t+1}|y_t,\ldots,y_1), \\ +F_{t} &= \VAR(v_t) = Z_t P_t Z_t^\top + H_t. +\end{aligned} +\end{equation*} + +Using the results of the Kalman filtering, we establish the state smoothing equations running backwards in time and yielding +\begin{equation*} +\begin{aligned} +\hat\alpha_{t} &= \E(\alpha_{t}|y_n,\ldots,y_1), \\ +V_{t} &= \VAR(\alpha_{t}|y_n,\ldots,y_1). +\end{aligned} +\end{equation*} + +Similar smoothed estimates can also be computed for the disturbance terms $\epsilon_t$ and $\eta_t$, and straightforwardly for the mean signal $\theta_t=Z_t\alpha_t$. For details on these algoritms, see Appendix \ref{appendix} and \citet{DK2012}. + +A prior distribution of the initial state vector $\alpha_1$ can be defined as a multivariate Gaussian distribution with mean $a_1$ and covariance matrix $P_1$. For uninformative diffuse prior, one typically sets $P_1 = \kappa\textrm{I}$, where $\kappa$ is for example $10^7$. However, this method can be numerically unstable due to cumulative roundoff errors. To solve this issue \citet{KD2003} present the exact diffuse initialization method, where the diffuse elements in $a_1$ are set to zero and $P_1$ is decomposed as $\kappa P_{\infty,1} + P_{*,1}$, where $\kappa\to\infty$. Here $P_{\infty,1}$ is a diagonal matrix with ones on those diagonal elements which relate to the nonstationary elements of $\alpha_1$, and $P_{*,1}$ contains the covariances of the stationary elements of $\alpha_1$ (and zeros elsewhere). At the start of the Kalman filtering (and at the end of backward smoothing) we use so called exact diffuse initialisation formulas until $P_{\infty,t}$ becomes zero matrix, and then continue with usual Kalman filtering equations. This exact method should be less prone to numerical errors, although they can still occur especially in the smoothing phase, if we have for example high collinearity between the explanatory variables of the model. Note that given all the parameters in the system matrices, results from the Kalman filter and smoother are equivalent with Bayesian analysis given the same prior distribution for $\alpha_1$. + +When we have multivariate observations, it is possible that in the diffuse phase, the matrix $F_t$ is not +invertible, and the computation of $a_{t+1}$ and $P_{t+1}$ becomes +impossible. On the other hand, even if $F_t$ is invertible, the computations can +become slow when dimensionality of $F_t$, i.e., the number of series increases. Also in case of multivariate observations, the formulas relating to the diffuse initialization become cumbersome. Based on the ideas of \citet{AM1979}, a complete univariate approach for filtering and smoothing was introduced by \citet{KD2000} (known as sequential processing by Anderson and Moore). The univariate approach is based on the alternative representation of the model \eqref{ssgeneral}, namely +\begin{equation*} +\begin{aligned} +y_{t,i} &= Z_{t,i}\alpha_{t,i} + \epsilon_{t,i}, \quad i=1,\ldots,p_t, \quad t=1,\ldots,n, \\ +\alpha_{t,i+1} &= \alpha_{t,i}, \quad i=1,\ldots,p_t-1,\\ +\alpha_{t+1,1} &= T_t\alpha_{t,p_t} + R_t\eta_{t}, \quad t=1,\ldots,n, +\end{aligned} +\end{equation*} +and $a_{1,1} \sim N(a_1,P_1)$, with an assumption that $H_t$ is diagonal for all $t$. Here the dimension of the observation vector $y_t$ can vary over time and therefore missing observations are handled +straightforwardly by adjusting the dimensionality of $y_t$. In case of +non-diagonal $H_t$, the original model can be transformed either by taking the +LDL decomposition of $H_t$, and multiplying the observation equation with the +$L_t^{-1}$, so $\epsilon_t^* \sim N(0,D_t)$, or by augmenting the state vector +with $\epsilon$, when $Q_t$ becomes block diagonal with blocks $Q_t$ and $H_t$. The augmenting can also be used for introducing correlation between $\epsilon$ and $\eta$. Both LDL decomposition and state vector augmentation is supported in \pkg{KFAS}. + +In theory when using the univariate approach, the computational costs of filtering and smoothing +decrease, as the number of matrix multiplications decrease, and there is no need +for solving the system of equations \citep[p. 159]{DK2012}. As noted in \citep{Tusell2010}, these gains can somewhat cancel out as more calls to linear algebra functions are needed and the memory management might not be as effective as working with larger objects at once. Nevertheless as noted previously, the sequential processing has also other clear benefits especially with diffuse initialization where the univariate approach simplifies the recursions considerably \citep{DK2012}. + +\pkg{KFAS} uses this univariate approach in all cases. Although $K_t=P_t Z_t^\top=\COV(a_t,y_t|y_t-1,\ldots,y_1)$, $v_t$, and $F_t$ differ from the standard multivariate versions, we get $a_t=a_{t,1}$ and $P_t=P_{t,1}$ by using the univariate approach. If standard multivariate matrices $F_t$ and $K_t$ are needed for inference, they can be computed later from the results of the univariate filter. As $F_{\ast,i,t}$, $K_{\ast,i,t}$, and $P_{\ast,t}$ coincide with the +nondiffuse counterparts if $F_{\infty,i,t}=0$, the asterisk is dropped from the variable names in \pkg{KFAS}, and for example variable \code{F} is a $n \times p$ array containing $F_{\ast,i,t}$ and $F_{i,t}$, whereas \code{Finf} is a $n \times d$, where $d$ is the last time point before the diffuse phase ended. + + +\subsection{Log-likelihood of the Gaussian state space model} + +The Kalman filter equations can be used for computing the log-likelihood, which +in its standard form is +\begin{equation*}\label{logL} +\log L = \frac{np}{2}\log2\pi - \frac{1}{2}\sum^n_{t=1}(\log|F_t| + +v^\prime_tF^{-1}_tv_t). +\end{equation*} +In case of the univariate treatment and diffuse initialization, the diffuse +log-likelihood can be written as +\begin{equation*}\label{logLd} +\begin{aligned} +\log L_d &= -\frac{1}{2}\sum^{n}_{t=1}\sum^{p_t}_{i=1}w_{i,t}, +\end{aligned} +\end{equation*} +where +\begin{displaymath} +w_{i,t} = \left\{ \begin{array}{ll} +\log F_{\infty,i,t}, &\mbox{if $F_{\infty,i,t}>0$,} \\ +I(F_{i,t}>0)(\log2\pi + \log F_{i,t} + v^2_{i,t}F_{i,t}^{-1}), &\mbox{if $F_{\infty,i,t}=0$}. +\end{array} \right. +\end{displaymath} + +See \citet[Chapter 7]{DK2012} for details. \citet{Francke2010} show that there are cases where the above definition of diffuse log-likelihood is not optimal. Without going into the details, if system matrices $Z_t$ or $T_t$ contain unknown parameters in their diffuse parts, the diffuse likelihood is missing one term which depends on those unknown parameters. \citet[p.411--412]{Francke2010} present a recursive formula for computing this extra term, which is also supported by \pkg{KFAS}. + +\section{State space models for exponential family}\label{exp} + +\pkg{KFAS} can also deal with observations which come from distributions of exponential family class other than Gaussian. We assume that the state equation as is in the Gaussian case, but the observation equation has the form +\begin{equation*} +p(y_t|\theta_t) = p(y_t|Z_t\alpha_t), +\end{equation*} +where $\theta_t=Z_t\alpha_t$ is the signal and $p(y_t|\theta_t)$ is the observational density. + +The signal $\theta_t$ is the linear predictor which is connected to the expected value $E(y_t)=\mu_t$ via a link function $l(\mu_t)=\theta_t$. In \pkg{KFAS}, the following distributions and links are available: + +\begin{enumerate} +\item Gaussian distribution with mean $\mu_t$ and variance $u_t$ with identity link $\theta_t=\mu_t$. + +\item Poisson distribution with intensity $\lambda_t$ and exposure $u_t$ together with log-link $\theta_t = \log(\lambda_t)$. Thus we have $\E(y_t|\theta_t)=\VAR(y_t|\theta_t)=u_t e^{\theta_t}$. + +\item Binomial distribution with size $u_t$ and probability of success $\pi_t$. \pkg{KFAS} uses logit-link so $\theta_t = \textrm{logit}(\pi_t)$ resulting $\E(y_t|\theta_t)=u_t\pi_t$ and $\VAR(y_t|\theta_t) = u_t(\pi_t(1-\pi_t))$. + +\item Gamma distribution with a shape parameter $u_t$ and an expected value $\mu_t$, again with log-link $\theta_t = \log(\mu)$, where Gamma distribution is defined as +\begin{equation*} +p(y_t|\mu_t,u_t) = \frac{u_t^{u_t}}{\Gamma(u_t)}\mu_t^{-u_t}y_t^{u_t-1}e^{\frac{y_t u_t}{\mu_t}}. +\end{equation*} +This gives us $\E(y_t|\theta_t)=e^{\theta_t}$ and $\VAR(y_t|\theta_t) = e^{2\theta_t}/u_t$. + +\item Negative binomial distribution with a dispersion parameter $u_t$ and an expected value $\mu_t$ with log-link $\theta_t = \log(\mu_t)$, where the negative binomial distribution is defined as +\begin{equation*} +p(y_t|\mu_t,u_t) = \frac{\Gamma(y_t+u_t)}{\Gamma(u_t)y_t!}\frac{\mu_t^{y_t} u_t^{u_t}}{(\mu_t+u_t)^{u_t+y_t}}, +\end{equation*} +giving us $\E(y_t|\theta_t)=e^{\theta_t}$ and $\VAR(y_t|\theta_t) = e^{\theta_t} + e^{2\theta_t}/u_t$. +\end{enumerate} + +Note that variable $u_t$ has a different meaning depending on the distribution it is linked to. In \pkg{KFAS} one defines the distribution for each time series via argument \code{distribution} and the additional known parameters $u_t$ corresponding to each series as columns of matrix \code{u}. + +In order to make inferences of the non-Gaussian models, we first find a Gaussian model which has the same conditional posterior mode as $p(\theta|y)$ \citep{DK2000}. This is done using an iterative process with Laplace approximation of $p(\theta|y)$, where the updated estimates for $\theta_t$ are computed via Kalman filtering and smoothing from the approximating Gaussian model. In approximating Gaussian model the observation equation is replaced by +\begin{equation*} +\tilde y_t = Z_t\alpha_t + \epsilon_t, \quad \epsilon_t \sim N(0,H_t) +\end{equation*} +where the pseudo-observations $\tilde y_t$ variances $H_t$ are based on first and second derivatives of $log p(y_t|\theta_t)$ with respect to $\theta_t$ \citep{DK2000}. + +Final estimates $\hat \theta_t$ correspond to the mode of $p(\theta|y)$. In Gaussian case the mode is also the mean. In cases listed in (1)-(5) the difference between the mode and mean is often negligible. Nevertheless, we are usually more interested in $\mu_t$ than in the linear predictor $\theta_t$. As the link function is non-linear, direct transformation $\hat \mu_t=l^{-1}(\hat \theta_t)$ introduces some bias. To solve this problem \pkg{KFAS} also contains methods based on importance sampling, which allows us to correct these possible approximation errors. With importance sampling technique we can also compute the log-likelihood and the smoothed estimates for $f(\alpha)$, where $f$ is an arbitrary function of states, $\exp(Z_t\alpha_t)$ being typical example. + +In importance sampling scheme, we first find the approximating Gaussian model, simulate the states $\alpha^i$ from this Gaussian model and then compute the corresponding weights $w_i=p(y|\alpha^i)/g(y|\alpha^i)$, where $p(y|\alpha^i)$ represents the conditional non-Gaussian density of the original observations, and $g(y|\alpha^i)$ is the conditional Gaussian density of the pseudo-observations $\tilde y$. These weights are then used for computing +\begin{equation*} +\E(f(\alpha)|y) = \frac{\sum_{i=1}^N f(\alpha^i) w_i}{\sum_{i=1}^Nw_i}. +\end{equation*} + +The simulation of Gaussian state space models in \pkg{KFAS} is based on simulation smoothing algorithm by \citet{DK2002}. In order to improve the simulation efficiency, \pkg{KFAS} can use two antithetic variables in the simulation algorithms. See \citet[p.~265-266]{DK2012} for details how these are constructed. + +\pkg{KFAS} also provides means for filtering of non-Gaussian models. This is achieved by sequentially using the smoothing scheme for $(y_1,\ldots,y_t), t=1\ldots,n$ with $y_t$ set as missing. This is relatively slow procedure for large models, as the importance sampling algorithms need to be performed $n$ times, although the first steps are much faster than the one using whole data. The non-Gaussian filtering is mainly for computation of recursive residuals (see Section \ref{residuals}) and for illustrative purposes, where the computational efficiency is not that important. With large models or online-filtering problems, one is recommended to use proper particle filter approach which is out of the scope of this paper. + +\subsection{Log-likelihood of the non-Gaussian state space model}\label{ngloglik} + +The log-likelihood function for the non-Gaussian model can be written as \citep[p.~272]{DK2012} +\begin{equation*} +\begin{aligned} +\log L(y) &= \log \int p(\alpha,y)\textrm{d}\alpha \\ + &= \log L_g(y)+ \log E_g\left[\frac{p(y|\theta)}{g(y|\theta)}\right], +\end{aligned} +\end{equation*} + +where $L_g(y)$ is the log-likelihood of the Gaussian approximating model and the expectation is taken with respect to the Gaussian density $g(\alpha|y)$. The expectation can be approximated by +\begin{equation}\label{eg} +\log E_g\left[\frac{p(y|\theta)}{g(y|\theta)}\right] \approx \log\frac{1}{N}\sum_{i=1}^N w_i. +\end{equation} + +In many cases, good approximation of the log-likelihood can be computed without any simulation, by setting $N=0$ and using the mode estimate $\hat \theta$ from the approximating model. + +In practice \eqref{eg} suffer from the fact that $w_i= p(y|\theta^i)/g(y|\theta^i)$ is numerically unstable; when number of observations is large, the discrete probability mass function $p(y|\theta^i)$ tends to zero, even when the Gaussian density function $g(y|\alpha^i)$ does not. Therefore it is better to redefine the weights as +\begin{equation*} +w^*_i = \frac{p(y|\theta^i)/p(y|\hat\theta)}{g(y|\theta^i)/g(y|\hat \theta)}. +\end{equation*} + +The log-likelihood is then computed as + +\begin{equation*} +\begin{aligned} +\log \hat L(y) &= \log L_g(y) + \log \hat w + \log\frac{1}{N}\sum_{i=1}^Nw^*_i, +\end{aligned} +\end{equation*} +where $\hat w = p(y|\hat \theta)/g(y|\hat \theta)$. + +\section{Residuals}\label{residuals} + +For exponential family state space models, multiple types of residuals can be computed. Probably the most useful ones are standardized recursive residuals, which are based on the one-step ahead predictions from the Kalman filter. For univariate case these are defined as +\begin{equation*} +\frac{y_t-\E(y_t|y_{t-1},\ldots,y_1)}{\sqrt{\VAR(y_t|y_{t-1},\ldots,y_1)}},\quad t = d+1\ldots,n, +\end{equation*} +where $d$ is the last time point of diffuse phase, and the denominator can be decomposed as +\begin{equation*} +\begin{aligned} +\VAR(y_t|y_{t-1},\ldots,y_1)&= \VAR(\E(y_t|\theta_t,y_{t-1},\ldots,y_1)|y_{t-1},\ldots,y_1)\\ +&+ \E(\VAR(y_t|\theta_t,y_{t-1},\ldots,y_1)|y_{t-1},\ldots,y_1)\\ +&= \VAR(\E(y_t|\theta_t)| y_{t-1},\ldots,y_1) + \E(\VAR(y_t|\theta_t)|y_{t-1},\ldots,y_1) . +\end{aligned} +\end{equation*} +In the Gaussian case this simplifies to $v_{t}F_{t}^{-\frac{1}{2}}$. + +For multivariate observations we have several options on how to standardize the residuals. The most common one is a marginal standardization approach where each residual series is divided by its standard deviation, so we get standard normal distributed residual series which have no autocorrelation or cross-correlation with non-zero lags. Another option is to use for example Cholesky decomposition for the prediction error covariance matrix $F_t$ and standardize the residuals by $L_t^{-1}(y_t-\hat y_t)$ where $L_t L_t^\top=F_t$. Now the whole residual series should look like a draw from standard normal distributed without any autocorrelation. + +For computing the marginally standardized residuals, multivariate versions of $F_t$ and $v_t$ are needed, wheras the Cholesky standardized residuals can be computed directly from sequential Kalman filter as +\begin{equation*} +v_{i,t} F_{i_t}^{-\frac{1}{2}}, \quad j=1,\ldots,p,\quad t=d+1\ldots,n. +\end{equation*} +These multivariate residuals depend on the ordering of the series, so if the residual diagnostics exhibit deviations from model assumptions, then the intepretation is somewhat more difficult than when using the marginal residuals. Therefore marginal residuals might be preferred. Note that if we want quadratic form residuals $(y_t- \hat y_t)F_t^{-1}(y_t- \hat y_t)$, then the ordering of the series does not matter. + +The recursive residuals are defined just for the non-diffuse phase, which is problematic if the model contains long diffuse phase for example if a dummy variable with a diffuse prior is incorporated to the model. This is because the diffuse phase cannot end before the dummy variable changes its value at least once. In order to circumvent this, one can set the a proper but highly non-informative prior distribution for the intervention variable when computing the residuals, which should have negligible effect on the visual inspection of the residual plots. + +Other potentially useful residuals are auxiliary residuals which are based on smoothed values of states. For details, see \citet{HarveyKoopman1992} and \citet[Chapter 7]{DK2012}. + +\section[Functionality of KFAS]{Functionality of \pkg{KFAS}} + +The state space model used with \pkg{KFAS} is built using function \code{SSModel}. The function uses \proglang{R}'s formula object in a similar way as for example functions \code{lm} and \code{glm}. In order to define the different components of the state space model, auxiliary functions \code{SSMtrend}, \code{SSMseasonal}, \code{SSMcycle}, \code{SSMarima}, \code{SSMregression} are provided These functions can be used to define the structural, ARIMA, and regression components of the model. The function \code{SSMcustom} can be used for constructing an arbitrary component by directly defining the system matrices of model \eqref{ssgeneral}. More details on how to construct common state space models with \pkg{KFAS} are presented in Section \ref{models}. + +The function \code{SSModel} returns an object of class \code{SSModel}, which contains the observations \code{y} as the \code{ts} object, system matrices \code{Z},\code{H},\code{T},\code{R},\code{Q} as arrays of appropriate dimensions, together with matrices \code{a1}, \code{P1}, and \code{P1inf} defining the initial state distribution. Additional components contains the system matrix \code{u} which is used in non-Gaussian models for additional parameters, character vector \code{distribution} which defines the distributions of the observations (multivariate series can have different distributions), and tolerance parameter \code{tol} which is used in diffuse phase for checking whether $F_\infty$ is nonzero. + +\code{SSModel} object also contains some attributes, namely integer valued attributes \code{p,m,k,n} which define the dimensions of the system matrices, and character vectors \code{state_types} and \code{eta_types} which define the elements of $\alpha_t$ and $\eta_t$. These attributes are used internally by \pkg{KFAS}, although user can carefully modify them if needed. For example, if the user wishes to redefine the error term $\eta_t$ by changing the dimensions of \code{R} and \code{Q}, the attributes \code{k} and \code{eta_types} need to be updated accordingly. + +The unknown model parameters can be estimated with \code{fitSSM}, which is a wrapper around the \pkg{R}'s \code{optim} function and the \code{logLik} method for \pkg{SSModel} object. For \code{fitSSM}, user gives the model object, initial values of unknown parameters and a function \code{updatefn} which is used to update the model given the parameters (the help page of \code{fitSSM} gives an example of \code{updatefn}). As the numerical optimization routines update the model and compute the likelihood thousands of times, the user is encouraged to build his own problem-specific model updating function for maximum efficiency. By default, \code{fitSSM} estimates the \code{NA} values in the time invariant covariance matrices $H$ and $Q$, but no general estimation function is provided. Of course, user can also directly use \code{logLik} method for computing the likelihood and thus is free to choose a suitable optimization method for his problem. + +Function \code{KFS} computes the filtered (one step ahead prediction) and smoothed estimates for states, signals, and the values of the inverse link function (expected value $\mu$ or probability $\pi$) in a non-Gaussian case. For Gaussian models, disturbance smoothing is also available. + +With \code{simulateSSM} user can simulate the states, signals or disturbances of the Gaussian state space models given the model and the observations. If the model contains missing observations, these can also be simulated by \code{simulateSSM} in similar way. It is also possible to simulate states from predictive distributions $p(\alpha_t|y_1,\ldots,y_{t-1})$, $t=1,\ldots,n$. For these simulations, instead of using marginal distributions $N(a_t,P_t)$, \pkg{KFAS} uses a modification of \citet{DK2002}, where smoothing is replaced by filtering. + +For non-Gaussian models, \code{importanceSSM} returns the states or signals simulated from the approximating Gaussian model, and the corresponding weights $w_i$, which can then be used to compute arbitrary functions of the states or signals. + +There are several \code{S3} methods available for \code{SSModel} and \code{KFS} objects. For both objects, simple \code{print} methods are provided, and for \code{SSModel} objects there is the \code{logLik} method. The \code{predict} method is for computation of the point predictions together with confidence or prediction intervals. Extraction operator \code{[} for extracting and replacing the subsets of model elements is also provided for class \code{SSModel}. Use of this method when modifying model is suggested instead of common list extractor \code{$}, as the latter can accidentally modify the dimensions of the corresponding model matrices. + +For \code{KFS} object, the methods \code{residuals}, \code{rstandard} and \code{hatvalues} are provided. Also, function \code{signal} can be used for extracting subsets of signals from \code{KFS} objects, for example the part of $Z_t\alpha_t$ that corresponds to the regression part of the model. + + +\section{Constructing common state space models with KFAS}\label{models} + +In this section we present some typical models which can be presented in a state space form. More examples can be found on the main help page of \pkg{KFAS} by typing \code{?KFAS} after the package is loaded via \code{library(KFAS)}. + +All the auxialiary functions used in formula argument of the function \code{SSModel} have some common arguments which are not directly related to the system matrices of the corresponding component. In complex multivariate models, important argument is \code{index}, which defines the series for which the corresponding component is constructed. For example, if we have four time series ($p=4$), we may want to use certain regression component only for series 2 and 4. In this case we use argument \code{index=c(2,4)} when calling the appropriate \code{SSMregression} function. By default the index is \code{1:p} so the component is constructed for all series. + +Another argument used in several auxiliary functions is \code{type}, which can take two possible values. Value \code{"distinct"} defines the component separately for each series defined by \code{index} (with covariance structure defined by argument \code{Q}), whereas value \code{"common"} constructs single component which applies to all series defined by \code{index}. For example we can define distinct local level components for all series together with covariance matrix which captures the dependencies of the different series, or we can define just a single local level component which is common to all series. + + +\subsection{Structural time series} + +Structural time series refers to class of state space models where the observed time series is decomposed into several underlying components, such as trend and seasonal effects. The basic structural time series model is of the form + +\begin{equation} +\begin{aligned}\label{structural} +y_t &= \mu_t + \gamma_t + c_t + \epsilon_t, \quad \epsilon_t \sim N(0, H_t),\\ +\mu_{t+1} &= \mu_t + \nu_t + \xi_t, \quad \xi_t \sim N(0, Q_{\textrm{level},t}), \\ +\nu_{t+1} &= \nu_t + \zeta_t, \quad \zeta_t \sim N(0, Q_{\textrm{slope},t}), +\end{aligned} +\end{equation} +where $\mu_t$ is the trend component, $\gamma_t$ is the seasonal component and $c_t$ is the cycle component. The seasonal component with period $s$ can be defined in a dummy variable form +\begin{equation*} +\gamma_{t+1} = -\sum_{j=1}^{s-1}\gamma_{t+1-j} + \omega_t, \quad \omega_t \sim N(0,Q_{\textrm{seasonal},t}), +\end{equation*} +or trigonometric form where +\begin{equation*} +\begin{aligned} +\gamma_{t} &= \sum_{j=1}^{\lfloor s/2 \rfloor}\gamma_{j,t}, \\ +\gamma_{j,t+1} &= \gamma_{j,t} \cos\lambda_j + \gamma^{\ast}_{j,t} \sin\lambda_j + \omega_{j,t}, \\ +\gamma^{\ast}_{j,t+1} &= - \gamma_{j,t} \sin\lambda_j + \gamma^{\ast}_{j,t} \cos\lambda_j + \omega^{\ast}_{j,t}, \quad j=1,\ldots, \lfloor s/2 \rfloor, +\end{aligned} +\end{equation*} +with $\omega_{j,t}$ and $\omega^{\ast}_{j,t}$ being independently distributed variables with $N(0, Q_{\textrm{seasonal},t})$ distribution and $\lambda_j = 2\pi j/s$. + +Cycle component with period $s$ is defined as + +\begin{equation*} +\begin{aligned} +c_{t+1} &= c_t\cos\lambda_c + c_t^\ast\sin\lambda_c + \omega_t,\\ +c^\ast_{t+1} &= - c_t\sin\lambda_c + c_t^\ast\cos\lambda_c + \omega^\ast_t, +\end{aligned} +\end{equation*} + +with $\omega_t$ and $\omega^{\ast}_{t}$ being independent variables from $N(0, Q_{\textrm{cycle},t})$ distribution and frequency $\lambda_c = 2\pi/s$. + +For non-Gaussian models the observation equation of \eqref{structural} is replaced by $p(y_t|\theta_t)$ where $\theta_t = \mu_t + \gamma_t + c_t$. Additional Gaussian noise term $\epsilon_t$ can also be included in $\theta_t$ using \code{SSMcustom} function (this is illustrated in Section \ref{illustration}) + +\code{SSModel} contains three auxiliary functions, \code{SSMtrend}, \code{SSMcycle}, and \code{SSMseasonal}, for building structural time series. Argument \code{degree} of \code{SSMtrend} defines the degree of the polynomial component, where \code{1} corresponds to local level model and \code{2} to local linear trend model. Higher order polynomials can also be defined with larger values. Another important argument for \code{SSMtrend} is \code{Q} which defines the covariance structure of the trend component. This is typically a list of $p \times p$ matrices (with $p$ being the number of series for which the component is defined), where the first matrix corresponds to level component ($\mu$ in \eqref{structural}), second to slope component $\nu$ and so forth. + +Function \code{SSMcycle} differs from \code{SSMtrend} only by one argument. \code{SSMcycle} does not have argument \code{degree}, but instead it has argument \code{period} which defines the length of the cycle $c_t$. Same argument is also used in function \code{SSMseasonal}, which contains also another important argument \code{sea.type}, which can be used to define whether user wants a \code{dummy} or \code{trigonometric} seasonal. + + + + +\subsection{ARIMA models} + +ARIMA models are another typical time series modelling framework, which are also possible to define as a state space model. Auxiliary \code{SSMarima} defines ARIMA model using vectors \code{ar} and \code{ma}, which define the autoregressive and moving average coefficients respectively. Function assumes that all series defined by the \code{index} have the same coefficients. Argument \code{d} defines the degree of differencing, and logical argument \code{stationary} defines whether stationarity (after differencing) is assumed (if not, a diffuse initial states are used instead of stationary distribution). Univariate ARIMA($p$,$d$,$q$) model can be written as +\begin{equation*} +y^*_t = \phi_1 y^*_{t-1} + \ldots + \phi_p y^*_{t-p} + \xi_t + \theta_1 \xi_{t-1} + \ldots + \theta_q \xi_{t-q}, +\end{equation*} +where $y^*_t = \Delta^d y_t$ and $\xi_t \sim \textrm{N}(0,\sigma^2)$. Let $r=\max(p,q+1)$. \pkg{KFAS} defines the state space representation of ARIMA($p$,$d$,$q$) model with stationary initial distribution as +\begin{displaymath} +Z_t^\top = +\left( \begin{array}{c} +1_{d+1} \\ +0 \\ +\vdots \\ +0 +\end{array}\right), +T = +\left( \begin{array}{ccccc} +U_d & 1^\top_d & 0 & \cdots & 0\\ +0 &\phi_1 & 1 & & 0 \\ +\vdots & & & \ddots & \\ +\vdots & \phi_{r-1} & 0 & & 1 \\ +0 & \phi_{r} & 0 & \cdots & 0 +\end{array} \right), +R = +\left( \begin{array}{c} +0_{d} \\ +1 \\ +\theta_1 \\ +\vdots \\ +\theta_{r-1} +\end{array}\right), +\end{displaymath} +\begin{displaymath} +\alpha_t = +\left( \begin{array}{c} +y_{t-1} \\ +\vdots \\ +\Delta^{d-1} y_{t-1}\\ +y^*_t \\ +\phi_2 y^*_{t-1} + \ldots + \phi_r y^*_{t-r+1} + \theta_1\eta_t+\ldots+\theta_{r-1}\eta_{t-r+2} \\ +\vdots \\ +\phi_r y^*_{t-1} + \theta_{r-1}\eta_t +\end{array}\right), +Q=\sigma^2, +\end{displaymath} +\begin{displaymath} +a_1 = +\left( \begin{array}{c} +0 \\ +\vdots \\ +0 +\end{array}\right), +P_{*,1} = +\left( \begin{array}{cccccc} +0 & 0 \\ +0 & S_r +\end{array} \right), +P_{\infty,1} = +\left( \begin{array}{cccccc} +I_{d} & 0 \\ +0 & 0 +\end{array} \right), +\eta_t=\xi_{t+1} +\end{displaymath} +where $\phi_{p+1}=\ldots=\phi_r=\theta_{q+1}=\ldots=\theta_{r-1}=0$, $1_{d+1}$ is a $1\times (d+1)$ vector of ones, $U_d$ is $d \times d$ upper triangular matrix of ones and $S_r$ is the covariance matrix of stationary elements of $\alpha_1$. The elements of the initial state vector $\alpha_1$, which correspond to the differenced values $y_0, \ldots, \Delta^{d-1}y_0$ are treated as diffuse. The covariance matrix $S_r$ can be computed by solving the linear equation $(I-T \otimes T)\textrm{vec}(S_r) = \textrm{vec}(RR^\top)$ \citep[p.138]{DK2012}. + +Note that the \code{arima} function from \pkg{stats} also uses the same state space approach to univariate ARIMA models for estimating model coefficients. + + +\subsection{Linear and generalized linear models} + +An ordinary linear regression model +\begin{equation*} +y_t = x^\top_t \beta + \epsilon_t, \quad t=1,\ldots,n, +\end{equation*} +where $\epsilon_t \sim N(0,\sigma^2)$ can be written as a Gaussian state space model by defining $Z_t=x^\top_t$, $H_t=\sigma^2$, $R_t=Q_t=0$ and $\alpha_t=\beta$. Assuming that the prior distribution of $\beta$ is defined as diffuse, the diffuse likelihood of this state space model corresponds to a restricted maximum likelihood (REML). Then the estimate for $\sigma^2$ obtained from \code{fitSSM} would be the familiar unbiased REML estimate of residual variance. It is important to notice that for this simple model numerical optimization is not needed, since we can estimate $\sigma^2$ by running the Kalman filter with $H_t=1$, which gives us +\begin{equation*} +\hat \sigma^2 = \frac{1}{\sum I(F_{\infty,t} = 0)} \sum_{t=1}^n I(F_{\infty,t} = 0) v_t^2/F_t, +\end{equation*} +which equals to the REML estimate of $\sigma^2$. The initial Kalman filter already provides correct estimates of $\beta$ as $a_{n+1}$, and running Kalman filter again with $H_t=\sigma^2$ gives also the covariance matrix of $\hat \beta$ as $P_{n+1}$. + +The extension from linear model to generalized linear model is straightforward as the basic theory behind the exponential family state space modelling can be formulated from the theory of generalized linear models (GLM), and can be thought of as a extension to GLMs with additional dynamic structure. The iterative process of finding the approximating Gaussian model is equivalent with the famous iterative reweighted least squares (IRLS) algoritm \citep[p.~40]{MN1989}. If the model is ordinary GLM the final estimates of regression coefficients $\beta$ and their standard errors coincide with maximum likelihood estimates obtained from ordinary GLM fitting. By adjusting the prior distribution for $\beta$ we can use \pkg{KFAS} also for Bayesian analysis of Poisson and binomial regression (as those distributions do not depend on any additional parameters such as residual variance) with Gaussian prior. + +A simple (generalized) linear model can be defined using \code{SSModel} without any auxiliary functions by defining the regression formula in the main part of the \code{formula}. For example the following code defines a Poisson GLM which is identical to the one found in help page of \code{glm}: + +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlstd{counts} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{18}\hlstd{,}\hlnum{17}\hlstd{,}\hlnum{15}\hlstd{,}\hlnum{20}\hlstd{,}\hlnum{10}\hlstd{,}\hlnum{20}\hlstd{,}\hlnum{25}\hlstd{,}\hlnum{13}\hlstd{,}\hlnum{12}\hlstd{)} +\hlstd{outcome} \hlkwb{<-} \hlkwd{gl}\hlstd{(}\hlnum{3}\hlstd{,}\hlnum{1}\hlstd{,}\hlnum{9}\hlstd{)} +\hlstd{treatment} \hlkwb{<-} \hlkwd{gl}\hlstd{(}\hlnum{3}\hlstd{,}\hlnum{3}\hlstd{)} +\hlstd{d.AD} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(treatment, outcome, counts)} +\hlstd{glmModel1} \hlkwb{<-} \hlkwd{SSModel}\hlstd{(counts} \hlopt{~} \hlstd{outcome} \hlopt{+} \hlstd{treatment,} + \hlkwc{data} \hlstd{= d.AD,} \hlkwc{distribution} \hlstd{=} \hlstr{"poisson"}\hlstd{)} +\end{alltt} +\end{kframe} +\end{knitrout} + +The previous model could also be defined using the auxiliary function \code{SSMregression}: + +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlstd{glmModel2} \hlkwb{<-} \hlkwd{SSModel}\hlstd{(counts} \hlopt{~} + \hlkwd{SSMregression}\hlstd{(}\hlopt{~}\hlstd{outcome} \hlopt{+} \hlstd{treatment,} \hlkwc{data} \hlstd{= d.AD),} + \hlkwc{distribution} \hlstd{=} \hlstr{"poisson"}\hlstd{)} +\end{alltt} +\end{kframe} +\end{knitrout} + +Note also the \code{data} argument in \code{SSMregression}. This overrides the \code{data} argument of \code{SSModel}, but both can exist at the same time. \code{SSModel} tries to be clever in finding the correct variables in case of multiple \code{data} arguments, see example in the help page of \code{SSModel} for an illustration. + +If our observations are multivariate, a distinct regression components are defined for each of the series. For example, if counts \code{counts} above were a bivariate series, then both series would have own regression coefficients but same covariate values. By using \code{SSMregression} explicitly, one could also define \code{type="common"} which would construct common regression coefficients for all series. + +With \code{SSMregression} one can also define more complex regression models. The first argument of \code{SSMregression}, \code{rformula} can used to provide single formula or a list of formulas, where each component of the list contains the appropriate formula to be used for the corresponding series ($i$th formula in the list is used for the $i$th series defined by argument \code{index}). When \code{rformula} is a list, the \code{data} argument of \code{SSMregression} can be a single data frame (or environment), or a list of such data objects. If \code{data} is a list, $i$th element of that list is used for $i$th formula, and if \code{data} is a single data frame or environment, same data is used for all formulas. + +\subsection{Generalized linear mixed models} + +Just like in GLM setting, it is also possible to write generalized linear mixed model (GLMM) as a state space model. The difference between fixed and random effects lies in the initial state distribution; fixed effects are initialized via diffuse prior whereas random effects have proper variance defined by elements of $P_1$. Both types of states are automatically estimated by the Kalman filter, given the covariance structure of the random effects (and the residual variance or other parameters related to distribution of observation equation). + +In practice, the mixed model formulation becomes quite cumbersome especially in hierarchical settings, but with large longitudinal settings it might still be useful to write mixed model as state space model, as it is then straightforward to add for example stochastic cycles or trends to the model. An example of defining the linear mixed model using the sleep study data from \pkg{lme4} \citep{lme4article,lme4} package proceeds as follows: + +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlkwd{suppressWarnings}\hlstd{(}\hlkwd{library}\hlstd{(}\hlstr{"lme4"}\hlstd{,}\hlkwc{quietly}\hlstd{=}\hlnum{TRUE}\hlstd{))} +\hlcom{# Split data by grouping variable} +\hlstd{Y} \hlkwb{<-} \hlkwd{split}\hlstd{(sleepstudy[}\hlstr{"Reaction"}\hlstd{], sleepstudy[}\hlstr{"Subject"}\hlstd{])} +\hlstd{p} \hlkwb{<-} \hlkwd{length}\hlstd{(Y)} \hlcom{# Number of series} +\hlstd{Y} \hlkwb{<-} \hlkwd{matrix}\hlstd{(}\hlkwd{unlist}\hlstd{(Y),} \hlkwc{ncol} \hlstd{= p,} + \hlkwc{dimnames} \hlstd{=} \hlkwd{list}\hlstd{(}\hlkwa{NULL}\hlstd{,} \hlkwd{paste}\hlstd{(}\hlstr{"Subject"}\hlstd{,} \hlkwd{names}\hlstd{(Y))))} +\hlstd{dataf} \hlkwb{<-} \hlkwd{split}\hlstd{(sleepstudy, sleepstudy[}\hlstr{"Subject"}\hlstd{])} + +\hlcom{# Assume that we know the covariance structure } +\hlcom{# of random effects and the residual variance} + +\hlstd{covRandom} \hlkwb{<-} \hlkwd{matrix}\hlstd{(}\hlkwd{c}\hlstd{(}\hlnum{625}\hlstd{,}\hlnum{36}\hlstd{,}\hlnum{36}\hlstd{,}\hlnum{625}\hlstd{),}\hlnum{2}\hlstd{,}\hlnum{2}\hlstd{)} +\hlstd{sigma2} \hlkwb{<-} \hlnum{650} + +\hlcom{# Common fixed effect part and distinct random effect parts for each "series"} +\hlcom{# Set P1inf = 0 so diffuse initialization is not used for random effects} +\hlstd{lmmModel} \hlkwb{<-} \hlkwd{SSModel}\hlstd{(Y} \hlopt{~ -}\hlnum{1} + \hlopt{+} \hlkwd{SSMregression}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlkwd{list}\hlstd{(}\hlopt{~} \hlstd{Days), p),} \hlkwc{type} \hlstd{=} \hlstr{"common"}\hlstd{,} + \hlkwc{data} \hlstd{= dataf,} \hlkwc{remove.intercept} \hlstd{=} \hlnum{FALSE}\hlstd{)} + \hlopt{+} \hlkwd{SSMregression}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlkwd{list}\hlstd{(}\hlopt{~} \hlstd{Days), p),} \hlkwc{P1inf} \hlstd{=} \hlnum{0}\hlstd{,} + \hlkwc{data} \hlstd{= dataf,} \hlkwc{remove.intercept} \hlstd{=} \hlnum{FALSE}\hlstd{),} + \hlkwc{H} \hlstd{=} \hlkwd{diag}\hlstd{(sigma2, p))} +\hlcom{# Set covariance structure of the random effects which are states 3 to 38} +\hlcom{# One could also use more common way lmmModel$P1[-(1:2),-(1:2)] <- ...} +\hlstd{lmmModel[}\hlstr{"P1"}\hlstd{,} \hlnum{2}\hlopt{+}\hlnum{1}\hlopt{:}\hlstd{(}\hlnum{2}\hlopt{*}\hlstd{p)]} \hlkwb{<-} + \hlkwd{as.matrix}\hlstd{(}\hlkwd{.bdiag}\hlstd{(}\hlkwd{replicate}\hlstd{(p, covRandom,} \hlkwc{simplify} \hlstd{=} \hlnum{FALSE}\hlstd{)))} +\end{alltt} +\end{kframe} +\end{knitrout} + + + +\section{Illustration}\label{illustration} + +I now illustrate the use of \pkg{KFAS} with an example case. Our data consists of alcohol related deaths in Finland for years 1969--2012, in age groups 30--39, 40--49, 50--59 and 60--69. We also have an offset term of yearly population sizes in corresponding age groups. The data is taken from \citet{stat1, stat2}. As an illustration, we use only observations until 2007, and make predictions for years 2008--2013. Figure \ref{fig:alcoholPlot1} shows the number of deaths per 100,000 persons for all age groups. + +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlkwd{data}\hlstd{(}\hlstr{"alcohol"}\hlstd{)} +\hlkwd{colnames}\hlstd{(alcohol)} +\end{alltt} +\begin{verbatim} +## [1] "death at age 30-39" "death at age 40-49" +## [3] "death at age 50-59" "death at age 60-69" +## [5] "population by age 30-39" "population by age 40-49" +## [7] "population by age 50-59" "population by age 60-69" +\end{verbatim} +\begin{alltt} +\hlkwd{ts.plot}\hlstd{(}\hlkwd{window}\hlstd{(alcohol[,}\hlnum{1}\hlopt{:}\hlnum{4}\hlstd{]}\hlopt{/}\hlstd{alcohol[,}\hlnum{5}\hlopt{:}\hlnum{8}\hlstd{],} \hlkwc{end} \hlstd{=} \hlnum{2007}\hlstd{),} \hlkwc{col} \hlstd{=} \hlnum{1}\hlopt{:}\hlnum{4}\hlstd{,} + \hlkwc{ylab} \hlstd{=} \hlstr{"Alcohol related deaths in Finland per 100,000 persons"}\hlstd{,} + \hlkwc{xlab} \hlstd{=} \hlstr{"Year"}\hlstd{)} +\hlkwd{legend}\hlstd{(}\hlstr{"topleft"}\hlstd{,}\hlkwc{col} \hlstd{=} \hlnum{1}\hlopt{:}\hlnum{4}\hlstd{,} \hlkwc{lty} \hlstd{=} \hlnum{1}\hlstd{,} + \hlkwc{legend} \hlstd{=} \hlkwd{colnames}\hlstd{(alcohol)[}\hlnum{1}\hlopt{:}\hlnum{4}\hlstd{])} +\end{alltt} +\end{kframe}\begin{figure}[!ht] + +\includegraphics[width=\linewidth]{figure/alcoholPlot1} \caption[Alcohol related deaths per 100,000 persons in Finland in 1969--2007]{Alcohol related deaths per 100,000 persons in Finland in 1969--2007.\label{fig:alcoholPlot1}} +\end{figure} + + +\end{knitrout} + +Natural distributional assumption for modelling counts is a Poisson distribution. Based on the time series plot, we can think of several candidates for capturing the time series aspects of the series. One could try for example an ARIMA model or a structural time series model such as local level or local linear trend. These are closely related, but I feel that the latter models are more easily interpretable so I will use structural time series approach. + +Here I choose a multivariate Poisson model + +\begin{equation*}\label{alcohol} +\begin{aligned} +p(y_t|\theta_t) &= Poisson(u_t e^{\theta_t}), \quad u_t = \textrm{population}_t,\\ +\theta_t &= \mu_t + \epsilon_t, \quad \epsilon_t \sim N(0, Q_{\textrm{noise}}),\\ +\mu_{t+1} &= \mu_t + \nu_t + \xi_t, \quad \xi_t \sim N(0, Q_{\textrm{level}}), \\ +\nu_{t+1} &= \nu_t, +\end{aligned} +\end{equation*} + +where $\mu_t$ is the random walk with drift component, $\nu_t$ is a constant slope and $\epsilon_t$ is an additional white noise component which captures the extra variation of the series. I make no restrictions for the covariance structures of the level or the noise component. Note that the random walk with drift is a special case of local linear trend model where the covariance structure of the slope term is set to zero. + + +We estimate the model parameters first without simulation, and then using those estimates as initial values run the estimation procedure again with importance sampling. In this case, the results obtained from the importance sampling step are practically identical with the ones obtained from the initial step. + +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlcom{# remove the last observations} +\hlstd{alcoholPred} \hlkwb{<-} \hlkwd{window}\hlstd{(alcohol,} \hlkwc{start} \hlstd{=} \hlnum{1969}\hlstd{,} \hlkwc{end} \hlstd{=} \hlnum{2007}\hlstd{)} + +\hlstd{model} \hlkwb{<-} \hlkwd{SSModel}\hlstd{(alcoholPred[,}\hlnum{1}\hlopt{:}\hlnum{4}\hlstd{]} \hlopt{~} + \hlkwd{SSMtrend}\hlstd{(}\hlnum{2}\hlstd{,} \hlkwc{Q} \hlstd{=} \hlkwd{list}\hlstd{(}\hlkwd{matrix}\hlstd{(}\hlnum{NA}\hlstd{,}\hlnum{4}\hlstd{,}\hlnum{4}\hlstd{),} \hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{4}\hlstd{,}\hlnum{4}\hlstd{)))} \hlopt{+} + \hlkwd{SSMcustom}\hlstd{(}\hlkwc{Z} \hlstd{=} \hlkwd{diag}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{4}\hlstd{),} \hlkwc{T} \hlstd{=} \hlkwd{diag}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{4}\hlstd{),} + \hlkwc{Q} \hlstd{=} \hlkwd{matrix}\hlstd{(}\hlnum{NA}\hlstd{,}\hlnum{4}\hlstd{,}\hlnum{4}\hlstd{),} \hlkwc{P1} \hlstd{=} \hlkwd{matrix}\hlstd{(}\hlnum{NA}\hlstd{,}\hlnum{4}\hlstd{,}\hlnum{4}\hlstd{)),} + \hlkwc{distribution} \hlstd{=} \hlstr{"poisson"}\hlstd{,} \hlkwc{u} \hlstd{= alcoholPred[,}\hlnum{5}\hlopt{:}\hlnum{8}\hlstd{])} + +\hlcom{# Model updating function for fitSSM} +\hlstd{updatefn} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{pars}\hlstd{,} \hlkwc{model}\hlstd{,} \hlkwc{...}\hlstd{)\{} + \hlstd{Q} \hlkwb{<-} \hlkwd{diag}\hlstd{(}\hlkwd{exp}\hlstd{(pars[}\hlnum{1}\hlopt{:}\hlnum{4}\hlstd{]))} + \hlstd{Q[}\hlkwd{upper.tri}\hlstd{(Q)]} \hlkwb{<-} \hlstd{pars[}\hlnum{5}\hlopt{:}\hlnum{10}\hlstd{]} + \hlstd{model[}\hlstr{"Q"}\hlstd{,}\hlkwc{etas}\hlstd{=}\hlstr{"level"}\hlstd{]} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(Q)} + \hlstd{Q} \hlkwb{<-} \hlkwd{diag}\hlstd{(}\hlkwd{exp}\hlstd{(pars[}\hlnum{11}\hlopt{:}\hlnum{14}\hlstd{]))} + \hlstd{Q[}\hlkwd{upper.tri}\hlstd{(Q)]} \hlkwb{<-} \hlstd{pars[}\hlnum{15}\hlopt{:}\hlnum{20}\hlstd{]} + \hlstd{model[}\hlstr{"Q"}\hlstd{,}\hlkwc{etas}\hlstd{=}\hlnum{9}\hlopt{:}\hlnum{12}\hlstd{]} \hlkwb{<-} \hlstd{model[}\hlstr{"P1"}\hlstd{,}\hlkwc{states}\hlstd{=}\hlnum{9}\hlopt{:}\hlnum{12}\hlstd{]} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(Q)} + \hlstd{model} +\hlstd{\}} +\hlcom{# Initial the covariance structure of the random walks and extra noise} +\hlcom{# theta = log(intensity) = log(y/u)} +\hlcom{# covariance matrices are parameterized via log-Cholesky in fitSSM} +\hlstd{init} \hlkwb{<-} \hlkwd{chol}\hlstd{(}\hlkwd{cov}\hlstd{(}\hlkwd{log}\hlstd{(alcoholPred[,}\hlnum{1}\hlopt{:}\hlnum{4}\hlstd{]}\hlopt{/}\hlstd{alcoholPred[,}\hlnum{5}\hlopt{:}\hlnum{8}\hlstd{]))}\hlopt{/}\hlnum{10}\hlstd{)} + +\hlstd{fitinit} \hlkwb{<-} \hlkwd{fitSSM}\hlstd{(model,} \hlkwc{updatefn} \hlstd{= updatefn,} + \hlkwc{inits} \hlstd{=} \hlkwd{rep}\hlstd{(}\hlkwd{c}\hlstd{(}\hlkwd{log}\hlstd{(}\hlkwd{diag}\hlstd{(init)), init[}\hlkwd{upper.tri}\hlstd{(init)]),}\hlnum{2}\hlstd{),} + \hlkwc{method} \hlstd{=} \hlstr{"BFGS"}\hlstd{)} +\hlcom{# Now with simulation} +\hlstd{fit}\hlkwb{<-}\hlkwd{fitSSM}\hlstd{(model,} \hlkwc{updatefn} \hlstd{= updatefn,} + \hlkwc{inits} \hlstd{= fitinit}\hlopt{$}\hlstd{optim.out}\hlopt{$}\hlstd{par,} \hlkwc{method} \hlstd{=} \hlstr{"BFGS"}\hlstd{,} \hlkwc{nsim} \hlstd{=} \hlnum{250}\hlstd{)} +\hlstd{varcor} \hlkwb{<-} \hlstd{fit}\hlopt{$}\hlstd{model[}\hlstr{"Q"}\hlstd{,} \hlkwc{etas} \hlstd{=} \hlstr{"level"}\hlstd{]} +\hlstd{varcor[}\hlkwd{upper.tri}\hlstd{(varcor)]} \hlkwb{<-} \hlkwd{cov2cor}\hlstd{(varcor)[}\hlkwd{upper.tri}\hlstd{(varcor)]} +\hlkwd{print}\hlstd{(varcor,}\hlkwc{digits}\hlstd{=}\hlnum{2}\hlstd{)} \hlcom{#local level component} +\end{alltt} +\begin{verbatim} +## [,1] [,2] [,3] [,4] +## [1,] 0.0074 0.66022 0.8062 0.856 +## [2,] 0.0028 0.00239 0.1654 0.711 +## [3,] 0.0040 0.00047 0.0034 0.755 +## [4,] 0.0033 0.00156 0.0020 0.002 +\end{verbatim} +\begin{alltt} +\hlstd{varcor} \hlkwb{<-} \hlstd{fit}\hlopt{$}\hlstd{model[}\hlstr{"Q"}\hlstd{,} \hlkwc{etas} \hlstd{=} \hlstr{"custom"}\hlstd{]} +\hlstd{varcor[}\hlkwd{upper.tri}\hlstd{(varcor)]} \hlkwb{<-} \hlkwd{cov2cor}\hlstd{(varcor)[}\hlkwd{upper.tri}\hlstd{(varcor)]} +\hlkwd{print}\hlstd{(varcor,}\hlkwc{digits}\hlstd{=}\hlnum{2}\hlstd{)} \hlcom{#local level component #extra noise component} +\end{alltt} +\begin{verbatim} +## [,1] [,2] [,3] [,4] +## [1,] 0.00537 0.73118 0.75627 8.0e-01 +## [2,] 0.00315 0.00346 0.99924 9.9e-01 +## [3,] 0.00295 0.00313 0.00283 1.0e+00 +## [4,] 0.00043 0.00043 0.00039 5.4e-05 +\end{verbatim} +\begin{alltt} +\hlopt{-}\hlstd{fitinit}\hlopt{$}\hlstd{optim.out}\hlopt{$}\hlstd{val} \hlcom{#log-likelihood without simulation} +\end{alltt} +\begin{verbatim} +## [1] -704.8 +\end{verbatim} +\begin{alltt} +\hlopt{-}\hlstd{fit}\hlopt{$}\hlstd{optim.out}\hlopt{$}\hlstd{val} \hlcom{#log-likelihood with simulation } +\end{alltt} +\begin{verbatim} +## [1] -704.8 +\end{verbatim} +\end{kframe} +\end{knitrout} + +Parameter estimation of state space model is often a difficult task, as the likelihood surface contains multiple maxima, thus making the optimization problem highly dependent on the initial values. Often the unknown parameters are related to the unobserved latent states such as the covariance matrix in this example, without much a priori knowledge. Therefore, it is challenging to guess good initial values especially in more complex settings and multiple initial value configurations possibly with several different type of optimization routines is recommended before one can be reasonably sure that proper optimum is found. Here we use the covariance matrix of the observed series as initial values for the covariance structures. + +Another issue in case of non-Gaussian models is the fact that the likelihood computation is based on iterative procedure which is stopped using some stopping criteria (such as relative change of log-likelihood), so the function actually contains some noise. This in turn affects the gradient computations in BFGS and can in theory give unreliable results. Using derivative free method like Nelder-Mead is therefore sometimes recommended. On the other hand BFGS is usually much faster than Nelder-Mead and thus I prefer to try first BFGS at least in preliminary analysis. + +Using function \code{KFS} we can compute the smoothed estimates of states: + +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlstd{out} \hlkwb{<-} \hlkwd{KFS}\hlstd{(fit}\hlopt{$}\hlstd{model,} \hlkwc{nsim} \hlstd{=} \hlnum{1000}\hlstd{)} +\hlstd{out} +\end{alltt} +\begin{verbatim} +## +## Smoothed values of states and standard errors at time n = 39: +## Estimate Std. Error +## level.death at age 30-39 2.855916 0.078437 +## slope.death at age 30-39 0.010714 0.013714 +## level.death at age 40-49 4.031312 0.042376 +## slope.death at age 40-49 0.023719 0.007632 +## level.death at age 50-59 4.757803 0.039830 +## slope.death at age 50-59 0.050371 0.009585 +## level.death at age 60-69 4.493837 0.033290 +## slope.death at age 60-69 0.048239 0.007209 +## custom1 -0.000402 0.060395 +## custom2 -0.019549 0.040885 +## custom3 -0.016949 0.037024 +## custom4 -0.002134 0.005143 +\end{verbatim} +\end{kframe} +\end{knitrout} + +From the output of \code{KFS} we see that the slope term is not significant in the first age group. For time-varying states we can easily plot the estimated level and noise components, which shows clear trends in three age groups and highly correlated additional variation in all groups: + +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlkwd{plot}\hlstd{(}\hlkwd{coef}\hlstd{(out,}\hlkwc{states}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"level"}\hlstd{,}\hlstr{"custom"}\hlstd{)),} + \hlkwc{main} \hlstd{=} \hlstr{"Smoothed states"}\hlstd{,} \hlkwc{yax.flip}\hlstd{=}\hlnum{TRUE}\hlstd{)} +\end{alltt} +\end{kframe}\begin{figure}[!ht] + +\includegraphics[width=\linewidth]{figure/states} \caption[Smoothed level and white noise components]{Smoothed level and white noise components.\label{fig:states}} +\end{figure} + + +\end{knitrout} + +Note the large drop in noise component which relates to possible outlier in 1973 of the mortality series. As an illustration of model diagnostics, we compute recursive residuals for our model and check whether there is autocorrelation left in the residuals (Figure \ref{fig:diagnostics1}). + +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlstd{res} \hlkwb{<-} \hlkwd{rstandard}\hlstd{(}\hlkwd{KFS}\hlstd{(fit}\hlopt{$}\hlstd{model,} \hlkwc{filtering} \hlstd{=} \hlstr{"mean"}\hlstd{,} + \hlkwc{smoothing} \hlstd{=} \hlstr{"none"}\hlstd{,} \hlkwc{nsim} \hlstd{=} \hlnum{1000}\hlstd{))} +\hlkwd{acf}\hlstd{(res,} \hlkwc{na.action} \hlstd{= na.pass)} +\end{alltt} +\end{kframe}\begin{figure}[!ht] + +\includegraphics[width=\linewidth]{figure/diagnostics1} \caption[Autocorrelations and cross-correlations of recursive residuals]{Autocorrelations and cross-correlations of recursive residuals.\label{fig:diagnostics1}} +\end{figure} + + +\end{knitrout} + +We see occasional lagged cross-correlation between the residuals, but overall we can be relatively satisfied with our model. + +We can now predict the intensity $e^{\theta_t}$ of alcohol related deaths per 100,000 persons for each age group for years 2008--2015 using our estimated model. As our model is time varying (\code{u} varies), we need to provide the model for the future observations via \code{newdata} argument. In this case we can use \code{SSMcustom} function and provide all the necessary system matrices as once, together with constant \code{u=1} (our signal $\theta$ is already scaled properly as the original $u_t$ was the population per 100,000 persons). + +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlstd{pred}\hlkwb{<-}\hlkwd{predict}\hlstd{(fit}\hlopt{$}\hlstd{model,} \hlkwc{newdata} \hlstd{=} + \hlkwd{SSModel}\hlstd{(}\hlkwd{ts}\hlstd{(}\hlkwd{matrix}\hlstd{(}\hlnum{NA}\hlstd{,}\hlnum{6}\hlstd{,}\hlnum{4}\hlstd{),} \hlkwc{start} \hlstd{=} \hlnum{2008}\hlstd{)} \hlopt{~ -}\hlnum{1} + \hlopt{+} \hlkwd{SSMcustom}\hlstd{(}\hlkwc{Z} \hlstd{= fit}\hlopt{$}\hlstd{model}\hlopt{$}\hlstd{Z,} \hlkwc{T} \hlstd{= fit}\hlopt{$}\hlstd{model}\hlopt{$}\hlstd{T,} + \hlkwc{R} \hlstd{= fit}\hlopt{$}\hlstd{model}\hlopt{$}\hlstd{R,} \hlkwc{Q} \hlstd{= fit}\hlopt{$}\hlstd{model}\hlopt{$}\hlstd{Q),} \hlkwc{u} \hlstd{=} \hlnum{1}\hlstd{,} + \hlkwc{distribution}\hlstd{=} \hlstr{"poisson"}\hlstd{),} + \hlkwc{interval} \hlstd{=} \hlstr{"confidence"}\hlstd{,} \hlkwc{nsim} \hlstd{=} \hlnum{10000}\hlstd{)} +\end{alltt} +\end{kframe} +\end{knitrout} +\begin{knitrout} +\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe} +\begin{alltt} +\hlstd{trend} \hlkwb{<-} \hlkwd{exp}\hlstd{(}\hlkwd{signal}\hlstd{(out,} \hlstr{"trend"}\hlstd{)}\hlopt{$}\hlstd{signal)} +\hlkwd{par}\hlstd{(}\hlkwc{mfrow} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{2}\hlstd{,}\hlnum{2}\hlstd{),} \hlkwc{mar} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{2}\hlstd{,}\hlnum{2}\hlstd{,}\hlnum{2}\hlstd{,}\hlnum{2}\hlstd{)} \hlopt{+} \hlnum{0.1}\hlstd{,} \hlkwc{oma} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{2}\hlstd{,}\hlnum{2}\hlstd{,}\hlnum{0}\hlstd{,}\hlnum{0}\hlstd{))} +\hlkwa{for}\hlstd{(i} \hlkwa{in} \hlnum{1}\hlopt{:}\hlnum{4}\hlstd{)} + \hlkwd{ts.plot}\hlstd{(alcohol[,i]}\hlopt{/}\hlstd{alcohol[,}\hlnum{4}\hlopt{+}\hlstd{i],} + \hlstd{trend[,i],} + \hlstd{pred[[i]],} + \hlkwc{col} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{,}\hlkwd{rep}\hlstd{(}\hlnum{3}\hlstd{,}\hlnum{3}\hlstd{)),} \hlkwc{xlab} \hlstd{=} \hlkwa{NULL}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlkwa{NULL}\hlstd{,} + \hlkwc{main} \hlstd{=} \hlkwd{colnames}\hlstd{(alcohol)[i])} +\hlkwd{mtext}\hlstd{(}\hlstr{"Number of alcohol related deaths per 100,000 persons in Finland"}\hlstd{,} + \hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)} +\hlkwd{mtext}\hlstd{(}\hlstr{"Year"}\hlstd{,}\hlkwc{side}\hlstd{=}\hlnum{1}\hlstd{,}\hlkwc{outer}\hlstd{=}\hlnum{TRUE}\hlstd{)} +\end{alltt} +\end{kframe}\begin{figure}[!ht] + +\includegraphics[width=\linewidth]{figure/predictplot} \caption[Observed number of alcohol related deaths per 100,000 persons in Finland (black), fitted values (red) and intensity predictions for years 2008--2012 together with 95\% prediction intervals (green)]{Observed number of alcohol related deaths per 100,000 persons in Finland (black), fitted values (red) and intensity predictions for years 2008--2012 together with 95\% prediction intervals (green).\label{fig:predictplot}} +\end{figure} + + +\end{knitrout} + +Figure \ref{fig:predictplot} shows the observed deaths, smoothed trends for 1969--2007, and intensity predictions for 2008--2012 together with 95\% prediction intervals for intensity. When we compare our predictions to true observations, we see that in reality the number of deaths slightly increased in the oldest age group (ages 60--69), whereas in other age they decreased substantially during the forecasting period. This is partly explained by the fact that during this period the total alcohol consumption decreased almost monotonically, which in turn might have been caused by the increase in taxation of alcohol in 2008, 2009 and 2012. + +% \section[Comparison to INLA]{Comparison to \pkg{INLA}}\label{comparison} +% +% I will now shortly compare \pkg{KFAS} and \pkg{INLA}, which is an \proglang{R} package for approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximation \citep{Rue}. +% +% As an illustration, we reanalyze the salmonella data analyzed by \citet{Margolin} which is available from \pkg{INLA}. The data consists of the number of revertant colonies of TA98 Salmonella with different doses of quinoline. We model the number of colonies as a Poisson GLMM with two explanatory variables and a random intercept term which tries to capture the overdispersion in the data. The codes for inference with \pkg{INLA} used here can be found from \url{http://www.r-inla.org/examples/volume-1/code-for-salm-example}. +% +% <<'salmonellaINLA'>>= +% suppressWarnings(library("INLA",quietly = TRUE)) +% data("Salm") +% mod.salm <- inla(y ~ log(dose+10) + dose + +% f(rand, model = "iid", param = c(0.001,0.001)), +% family = "poisson", data = Salm) +% ## improved estimation of the hyperparameters +% h.salm <- inla.hyperpar(mod.salm) +% @ +% +% There are two ways to define the random intercept component in \pkg{KFAS}. The first one uses the \code{SSMregression} function and constructs a factor with 18 levels (one for each case) with non-diffuse initial variance $\sigma^2$. This gives 18 identically distributed time-invariant states where each state correspond to random effect of one observation. Another option would be to use \code{SSMcustom} function and define just one time-varying state as \code{SSMcustom(Z = 1, T = 0, R = 1, Q = sigma2, a1 = 0, P1 = sigma2, P1inf = 0)}. Both approaches give identical results. However, for large data the former approach is less efficient as the number of states depends on the number of observations. Nevertheless, we use the former approach here for illustration. +% +% <<'salmonellaKFAS'>>= +% Salm$rand<-as.factor(Salm$rand) +% model <- SSModel(y ~ log(dose+10) + dose + +% SSMregression(~ -1 + rand, P1 = diag(NA,18), +% remove.intercept = FALSE), +% data = Salm, distribution = "poisson") +% +% updatefn<-function(pars,model,...){ +% diag(model["P1", states = 4:21]) <- exp(pars) +% model +% } +% +% fit <- fitSSM(model, updatefn = updatefn, inits = -3, +% method = "BFGS", nsim = 1000) +% @ +% +% <<'salmonellaResults'>>= +% out<-KFS(fit$model,nsim=1000) +% out +% h.salm$summary.fixed[,1:2] +% h.salm$summary.random$rand[,2:3] +% 1/h.salm$summary.hyper[1] +% fit$model["P1", states = 4] +% @ +% +% Although \pkg{INLA} uses Bayesian approach which takes account the parameter estimation uncertainty, the results from \pkg{INLA} and \pkg{KFAS} are practically same, even with such a small data. Kalman filtering with diffuse initialization still takes account the uncertainty of the estimation of regression coefficients, so the differences here are related to the estimation of hyperparameter $\sigma^2$, which is estimated as precision $1/\sigma^2$ by \pkg{INLA}. The estimate of $\sigma^2$ by \pkg{KFAS} is signif(fit$model["P1", states = 4],2) whereas \pkg{INLA} gives $\sigma^2=signif(1/h.salm$summary.hyper[1],2)$. +% +% As \pkg{INLA} and \pkg{KFAS} are based on different (although related) theoretical framework, the extensive study of their performances in terms of computational efficiency and accuracy of results is somewhat pointless. Nevertheless, some remarks can be made. I feel that the biggest advantage of \pkg{INLA} is the Bayesian framework which allows us to take account the parameter uncertainty in predictions and other inference. On the other hand, the computational burden related to the numerical integration over the hyperparameters can become infeasible as the number of hyperparameters increase. It is not uncommon to have a time series model with tens (or even hundreds) of parameters (such as multivariate structural time series or dynamic factor models). Of course these same models cause problems also to the maximum likelihood estimation, as noted in the previous section. Also the Bayesian approach eliminates the need for defining good initial values for the maximum likelihood estimation but the problem transforms into defining good priors for the same hyperparameters which again is a non-trivial task in practice. + +\section{Discussion} + +State space models offers tools for solving a large class of statistical problems. Here I introduced an R package \pkg{KFAS} for linear state space modelling where the observations are from an exponential family. With such a general framework, different aspects of the modelling need to be taken into account. Therefore the focus of the package has been to provide reliable and relatively fast tools for multiple inference problems, such as maximum likelihood estimation, filtering, smoothing and simulation. Compared to the early versions of \pkg{KFAS}, constructing a state space model with simple components is now possible without explicit definition of the system matrices by using the auxiliary functions and symbolic descriptions with the help of formula objects, which should greatly ease the use of the package. + +Currently all the time consuming parts of \pkg{KFAS} are written in \proglang{Fortran}, which makes it relatively fast, given the general nature of problems \pkg{KFAS} can handle. Still, converting the package to \proglang{C++} and \code{S4} classes with help of \pkg{Rcpp} \citep{RcppA,RcppB} could result potential improvements in terms of memory management, scalability and maintenance. + +\section*{Acknowledgments} + +The author wishes to thank Jukka Nyblom and Patricia Menendez for the valuable comments and suggestions regarding the paper and the package. Author is also grateful for the financial support from Emil Aaltonen's Foundation. + +\clearpage +\newpage +\appendix +\section{Appendix: Filtering and smoothing recursions}\label{appendix} + +The following formulas summarize the Kalman filtering and smoothing formulas for diffuse and sequential case and are based on \citet{DK2001} and related articles. The original formulas are somewhat scattered between the references with slightly different notations. Therefore I have collected the equations used in \pkg{KFAS} to this Appendix. +\subsection{Filtering} +Denote +\begin{equation*} +\begin{aligned} +a_{t+1} &= \E(\alpha_{t+1}|y_t,\ldots,y_1) \quad \textrm{and} \\ +P_{t+1} &= \VAR(\alpha_{t+1}|y_t,\ldots,y_1). +\end{aligned} +\end{equation*} +The Kalman filter recursions for the general Gaussian model of form \eqref{ssgeneral} are +\begin{equation*} +\begin{aligned} +v_t &= y_t - Z_t a_t \\ +F_t &= Z_t P_t Z_t^\top + H_t \\ +K_t &= P_t Z_t^\top \\ +a_{t+1} &= T_t (a_t + K_t F^{-1}_t v_t) \\ +P_{t+1} &= T_t (P_t -K_tF^{-1}_tK_t^\top)T_t^\top + R_t Q_t R_t, +\end{aligned} +\end{equation*} + +For the univariate approach, the filtering equations are +\begin{equation*} +\begin{aligned} +v_{t,i} &= y_{t,i} - Z_{t,i} a_{t,i} \\ +F_{t,i} &= Z_{t,i} P_{t,i} Z_{t,i}^\top + \sigma^2_{t,i} \\ +K_{t,i} &= P_{t,i} Z_{t,i}^\top \\ +a_{t,i+1} &= a_{t,i} + K_{t,i}F_{t,i}^{-1} v_{t,i} \\ +P_{t,i+1} &= P_{t,i} -K_{t,i}K^\top_{t,i}F_{t,i}^{-1}\\ +a_{t+1,1} &= T_t a_{t,p_t+1}\\ +P_{t+1,1} &= T_t P_{t,p_t+1} T_t^\top + R_t Q_t R_t, +\end{aligned} +\end{equation*} +for $t=1,\ldots,n$ and $i=1,\ldots,p_t$, where $v_{t,i}$ and $F_{t,i}$ are +scalars, $K_{t,i}$ is a column vector and $\sigma^2_{t,i}$ is the $i$th diagonal +element of $H_t$. It is possible that $F_{t,i}=0$, which case $a_{t,i+1} = +a_{t,i}$, $P_{t,i+1} = P_{t,i}$, and $v_{t,i}$ is computed as usual. + +The diffuse filtering equations for univariate approach are +\begin{equation*} +\begin{aligned} +v_{t,i} &= y_{t,i} - Z_{t,i} a_{t,i} \\ +F_{\ast,t,i} &= Z_{t,i} P_{\ast,t,i} Z_{t,i}^\top + \sigma^2_{t,i} \\ +F_{\infty,t,i} &= Z_{t,i} P_{\infty,t,i} Z_{t,i}^\top \\ +K_{\ast,t,i} &= P_{\ast,t,i} Z_{t,i}^\top \\ +K_{\infty,t,i} &= P_{\infty,t,i} Z_{t,i}^\top, +\end{aligned} +\end{equation*} + +and + +\begin{equation*} +\begin{aligned} +a_{t,i+1} &= a_{t,i} + K_{\infty,t,i} v_{t,i}F_{\infty,t,i}^{-1} \\ +P_{\ast,t,i+1} &= P_{\ast,t,i} ++K_{\infty,t,i}K^\top_{\infty,t,i}F_{\ast,t,i}F_{\infty,t,i}^{-2} +-(K_{\ast,t,i}K^\top_{\infty,t,i}+K_{\ast,t,i}K^\top_{\infty,t,i})F_{\infty,t,i}^{-1}\\ +P_{\infty,t,i+1} &= P_{\infty,t,i} +-K_{\infty,t,i}K^\top_{\infty,t,i}F_{\infty,t,i}^{-1}\\ +\end{aligned} +\end{equation*} + +if $F_{\infty,t,i}>0$, and + +\begin{equation*} +\begin{aligned} +a_{t,i+1} &= a_{t,i} + K_{\ast,t,i} v_{t,i}F_{\ast,t,i}^{-1} \\ +P_{\ast,t,i+1} &= P_{\ast,t,i} -K_{\ast,t,i}K^\top_{\ast,t,i}F_{\ast,t,i}^{-1}\\ +P_{\infty,t,i+1} &= P_{\infty,t,i},\\ +\end{aligned} +\end{equation*} + +if $F_{\infty,t,i}=0$. The transition equations from $t$ to $t+1$ are + +\begin{equation*} +\begin{aligned} +a_{t+1,1} &= T_t a_{t,p_t+1}\\ +P_{\ast,t+1,1} &= T_t P_{\ast,t,p_t+1} T_t^\top + R_t Q_t R_t\\ +P_{\infty,t+1,1} &= T_t P_{\infty,t,p_t+1} T_t^\top. +\end{aligned} +\end{equation*} + + +\subsection{Smoothing} +Denote +\begin{equation*} +\begin{aligned} +\hat\alpha_{t} &= \E(\alpha_{t}|y_n,\ldots,y_1) \quad \textrm{and} \\ +V_{t} &= \VAR(\alpha_{t}|y_n,\ldots,y_1). +\end{aligned} +\end{equation*} + +The smoothing algorithms of \pkg{KFAS} are based on the following recursions: + +\begin{equation*} +\begin{aligned} +r_{t,i-1} &= Z_{t,i}^\top v_{t,i}F_{t,i}^{-1} + L_{t,i}^\top r_{t,i}, \\ +r_{t-1,p_t} &= T_{t-1}^\top r_{t,0},\\ +N_{t,i-1} &= Z_{t,i}^\top Z_{t,i}F_{t,i}^{-1} +L_{t,i}^\top N_{t,i}L_{t,i},\\ +N_{t-1,p_t} &= T_{t-1}^\top N_{t,0}T_{t-1},\\ +L_{t,i} &= I - K_{t,i}Z_{t,i}^\top F_{t,i}^{-1}, +\end{aligned} +\end{equation*} + +for $t=n,\ldots,1$ and $i=p_t,\ldots,1$, with $r_{n,p_n}=0$ and $N_{n,p_n}=0$. +From these recursions, we get state smoothing recursions +\begin{equation*} +\begin{aligned} +\hat\alpha_{t} &= a_{t,1} + P_{t,1}r_{t,0}\\ +V_{t} &= P_{t,1}-P_{t,1}N_{t,0}P_{t,1}, +\end{aligned} +\end{equation*} + +and disturbance smoothing recursions + +\begin{equation*} +\begin{aligned} +\hat \epsilon_{t,i} &= \sigma^2_{t,i}F_{t,i}^{-1}(v_{t,i} -K_{t,i}^\top r_{t,i}),\\ +\VAR(\hat \epsilon_{t,i}) &= \sigma^2_{t,i} - +\sigma^4_{t,i}(F_{t,i}^{-1} -K_{t,i}^\top N_{t,i}K_{t,i}F_{t,i}^{-2}),\\ +\hat \eta_{t} &= Q_tR_t^\top r_{t,0},\\ +\VAR(\hat \eta_{t,i}) &= Q_tR^\top_t N_{t,0}R_tQ_t. +\end{aligned} +\end{equation*} +The recursions for diffuse phase are as follows. +\begin{equation*} +\begin{aligned} +L_{\infty,t,i} &= I - K_{\infty,t,i}Z_{t,i}F_{\infty,t,i}^{-1},\\ +L_{t,i} &= +(K_{\infty,t,i}F_{t,i}F_{\infty,t,i}^{-1}-K_{t,i})Z_{t,i}F_{\infty,t,i}^{-1},\\ +r_{0,t,i-1} &= L_{\infty,t,i}^\top r_{0,t,i}, \\ +r_{1,t,i-1} &= Z_{t,i}^\top v_{t,i}F_{\infty,t,i}^{-1} + L_{\infty,t,i}^\top r_{1,t,i} + +L^\top_{t,i}r_{0,t,i}, \\ +N_{0,t,i-1} &= L_{\infty,t,i}^\top N_{0,t,i}L_{\infty,t,i}\\ +N_{1,t,i-1} &= L_{t,i}^\top N_{0,t,i}L_{\infty,t,i}+ +L_{\infty,t,i}^\top N_{1,t,i}L_{\infty,t,i}+Z_{t,i}^\top Z_{t,i}F_{\infty,t,i}^{-1},\\ +N_{2,t,i-1} &= L_{t,i}^\top N_{0,t,i}L_{t,i}+ L_{\infty,t,i}^\top N_{1,t,i}L_{t,i} + +(L_{\infty,t,i}^\top N_{1,t,i}L_{t,i})^\top + +L_{\infty,t,i}N_{2,t,i}^\top L_{\infty,t,i}-Z_{t,i}^\top Z_{t,i}F_{t,i}F_{\infty,t,i}^{-2},\\ +N_{t-1,p_t} &= T_{t-1}^\top N_{t,0}T_{t-1}, +\end{aligned} +\end{equation*} +if $F_{\infty,t,i}>0$, and +\begin{equation*} +\begin{aligned} +L_{t,i} &= I - K_{t,i}Z_{t,i}F_{t,i}^{-1},\\ +r_{0,t,i-1} &= Z_{t,i}^\top v_{t,i}F_{t,i}^{-1} + L_{t,i}^\top r_{0,t,i}, \\ +r_{1,t,i-1} &= L_{t,i}^\top r_{1,t,i}, \\ +N_{0,t,i-1} &= L_{t,i}^\top N_{0,t,i}L_{t,i}+Z_{t,i}^\top Z_{t,i}F_{t,i}^{-1}\\ +N_{1,t,i-1} &= N_{1,t,i}L_{t,i}\\ +N_{2,t,i-1} &= N_{2,t,i}L_{t,i}, +\end{aligned} +\end{equation*} +otherwise. The transition from time $t$ to $t-1$ is by $N_{j,t-1,p_t} = +T_{t-1}^\top N_{j,t,0}T_{t-1}$ for $j=0,1,2$, and $r_{j,t-1,p_t} = T_{t-1}^\top r_{j,t,0}$ +for $j=0,1$, with $r_{0,d,j}=r_{d,j}$, $r_{1,d,j}=0$, $N_{0,d,j}=N_{d,j}$, and +$N_{1,d,j}=N_{2,d,j}=0$, where $(d,j)$ is the last point of diffuse phase. +From these basic recursions, we get state smoothing recursions for diffuse phase +as +\begin{equation*} +\begin{aligned} +\hat\alpha_{t} &= a_{t,1} + P_{t,1}r_{0,t,0} + P_{\infty,t,1}r_{1,t,0},\\ +V_{t} &= P_{t,1}-P_{t,1}N_{0,t,0}P_{t,1} -(P_{\infty,t,1}N_{1,t,0}P_{t,1})^\top +-P_{\infty,t,1}N_{1,t,0}P_{t,1} - P_{\infty,t,1}N_{2,t,0}P_{\infty,t,1}, +\end{aligned} +\end{equation*} +and disturbance smoothing recursions +\begin{equation*} +\begin{aligned} +\hat \epsilon_{t,i} &= -\sigma^2_{t,i}K_{\infty,t,i}^\top r_{0,t,i},\\ +\VAR(\hat \epsilon_{t,i}) &= \sigma^2_{t,i} - +\sigma^4_{t,i}K_{\infty,t,i}^\top N_{0,t,i}K_{\infty,t,i}F_{\infty,t,i}^{-2}, +\end{aligned} +\end{equation*} +if $F_{\infty,t,i}>0$, and + +\begin{equation*} +\begin{aligned} +\hat \epsilon_{t,i} &= +-\sigma^2_{t,i}(v_{t,i}F_{\infty,t,i}^{-1}-K_{t,i}^\top r_{0,t,i}),\\ +\VAR(\hat \epsilon_{t,i}) &= \sigma^2_{t,i} - +\sigma^4_{t,i}(F_{t,i}^{-1}-K_{t,i}^\top N_{0,t,i}K_{t,i}F_{t,i}^{-2}), +\end{aligned} +\end{equation*} + +if $F_{\infty,t,i}=0$. For $\hat\eta$, recursions are + +\begin{equation*} +\begin{aligned} +\hat \eta_{t} &= Q_tR_t^\top r_{0,t,0},\\ +\VAR(\hat \eta_{t,i}) &= Q_tR_t^\top N_{0,t,0}R_tQ_t. +\end{aligned} +\end{equation*} + + + +\bibliography{KFAS} +\end{document} + diff --git a/vignettes/figure/alcoholPlot1.pdf b/vignettes/figure/alcoholPlot1.pdf new file mode 100644 index 0000000..4127bb3 Binary files /dev/null and b/vignettes/figure/alcoholPlot1.pdf differ diff --git a/vignettes/figure/diagnostics1.pdf b/vignettes/figure/diagnostics1.pdf new file mode 100644 index 0000000..55977fd Binary files /dev/null and b/vignettes/figure/diagnostics1.pdf differ diff --git a/vignettes/figure/predictplot.pdf b/vignettes/figure/predictplot.pdf new file mode 100644 index 0000000..c40108b Binary files /dev/null and b/vignettes/figure/predictplot.pdf differ diff --git a/vignettes/figure/states.pdf b/vignettes/figure/states.pdf new file mode 100644 index 0000000..069a51e Binary files /dev/null and b/vignettes/figure/states.pdf differ