! ARMAFILTER
! Performs ARMA(P,Q) filter
! SYNTAX:
!
! e = ARMAfilter(rho,psi,mu,Y,Q,P)
!
! e(t) = Y(t) - mu(t) - rho(1)Y(t-1) -...- rho(P)Y(t-P) 
!          - psi(1)e(t-1) - ... - psi(Q)e(t-Q)
! 
! Accompanying the GARCHKIT, there should be a mex-file version of 
! this function, as well as the fortran code used to create it. 
! This m-code has been provided in case the user can't or won't utilize
! the compiled code.
! Copyright (C) 2000 Matthew C. Roberts, matthewcroberts@hotmail.com 

! This toolbox is free software; feel free to modify, fix or improve it.
! Please, however, upon making substantive changes that you
! feel improve the software, send me a copy. If I include your version,
! I will credit you. Note, however, that these functions are free. You
! may freely make use of them, but may not charge for their distribution.
! 
! Finally, if a friend or colleague would like to use GARCHKIT, 
! please refer them to me at the address posted above. This way,
! they will have the latest revision, and i will be able
! to make them aware of future updates.
!
! This library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 

! =============================================================
! GATEWAY ROUTINE:
! SYNTAX (FROM MATLAB:)
! e = ARMAfilter(rho,psi,mu,Y,Q,P)
! =============================================================

!============================================================================
!	-----------------------------------------------------------------------
!	This is the calling function
!	-----------------------------------------------------------------------
subroutine mexFunction(nlhs,plhs,nrhs,prhs)
implicit none
integer  nrhs, nlhs
integer  m(nrhs), n(nrhs), i, j

integer  K, P, Q, T
real*8, allocatable :: mu(:,:), psi(:,:), rho(:,:), e(:,:), &
                       Y(:), tmpe(:), tmpY(:)
real*8 :: tmpReal

!	THE FOLLOWING SHOULD BE CHANGED TO INTEGER*8 ON 64-BIT MACHINES.	
integer mxIsNumeric, mxGetM, mxGetN, prhs(nrhs), plhs(nrhs), & 
    mxCreateFull, x(nrhs), outpr, mxGetPr

external mxIsNumeric, mxGetPr, mxCreateFull, mxGetM, mxGetN

!	Check for proper number of calling arguments
if (nlhs .gt. 1) then
	call mexErrMsgTxt('Only one return argument')
elseif (nrhs .ne. 6) then
	call mexErrMsgTxt('Must have exactly 6 calling arguments')
endif

!	Make sure the matrix is numeric
do i=1,nrhs
	if (mxIsNumeric(prhs(i)) .eq. 0) then
		call mexErrMsgTxt('Non-numeric data present')
	endif
	m(i)=mxGetM(prhs(i))
	n(i)=mxGetN(prhs(i))
end do

!	Make sure that if the matrix if vectorized, it all
!	all conforms properly
do j=1,3
  if (n(j) .ne. n(1)) then
	call mexErrMsgTxt('Inputs must have same number of columns')
  endif
end do

!	Make sure that inputs 5 & 6 are scalars:
if ((m(5) .ne. 1) .or. (n(5) .ne. 1)) then
  call mexErrMsgTxt('P must be a scalar')
end if
if ((m(6) .ne. 1) .or. (n(6) .ne. 1)) then
  call mexErrMsgTxt('P must be a scalar')
end if

!	Get pointers:
do i=1,nrhs
	x(i) = mxGetPr(prhs(i))
end do

!	rename the size params:
K  = n(1)
P  = m(1)
Q  = m(2)
T  = m(4)

! allocate matrices to hold the input & output args:
allocate( mu(T,K), psi(Q,K), rho(P,K), e(T,K), Y(T) )
allocate( tmpe(T+max(P,Q)), tmpY(T+max(P,Q)) )

!	Copy data to local workspaces:
if (P*K .gt. 0) then
  call mxCopyPtrtoReal8(x(1),rho,P*K)
else
  rho = 0.0d0
end if
if (Q*n(2) .gt. 0) then
  call mxCopyPtrtoReal8(x(2),psi,Q*K)
else
  psi = 0.0d0
end if
call mxCopyPtrtoReal8(x(3),mu,T*K)
call mxCopyPtrtoReal8(x(4),Y,T)

!	Prepare for output function:
plhs(1)=mxCreateFull(T,K,0.0d0)
outpr=mxGetPr(plhs)

!     Y into tmpY:
tmpY(1+max(P,Q):) = Y

!     set e=0
tmpe = 0.0d0

!	call numerical function:
do j = 1,K
  tmpReal = sum(rho(:,j)) + sum(psi(:,j))
  tmpY(1:max(P,Q)) = mu(1,j)/(1.0d0 - tmpReal)
  call ARMAfilter(rho(:,j),psi(:,j),mu(:,j),tmpY,P,Q,T,tmpe,e(:,j))
end do
!	copy results out:
call mxCopyReal8toPtr(e,outpr,T*K)

contains

! ======================================================
!	Here is the numerical routine!!!!
!      e = ARMAfilter(rho,psi,mu,Y,Q,P)
! ======================================================
subroutine ARMAfilter(rho,psi,mu,tmpY,P,Q,T,tmpe,e)
implicit none
integer P, Q, Tm, T
double precision rho(P), psi(Q), mu(T), tmpe(T+max(P,Q)), &
    e(T), tmpY(T+max(P,Q))
 
integer i, j, MaxPQ

MaxPQ = max(P,Q)

do i=1,T
  tmpe(i+MaxPQ) = tmpY(i+MaxPQ) - mu(i)
  do j=1,P
	tmpe(i+MaxPQ) = tmpe(i+MaxPQ) - rho(j)*tmpY(i+MaxPQ-j)
  end do
  do j=1,Q
	tmpe(i+MaxPQ) = tmpe(i+MaxPQ) - psi(j)*tmpe(i+MaxPQ-j)
  end do
end do

e = tmpe(MaxPQ+1:)


end subroutine ARMAfilter

end subroutine MexFunction