module altmod	! likelihood under alternative where p=p1
	use likmod
	implicit none
	
	contains
	
	function getstatvar1(th1) result(val)
		use lincg_int
		use iercd_int
		real			:: th1(nth1),val(p1,p1)
		complex	:: lam(p1), Q(p1,p1), Qinv(p1,p1)
		complex	:: a(p1**2), b(p1**2)
		integer	:: i,j
		lam=-getcfromh(th1(1:p1))
		do j=1,p1
			Q(1,j)=1.0
			do i=2,p1
				Q(i,j)=lam(j)*Q(i-1,j)
			enddo
			Q(:,j)=Q(:,j)/norm2(abs(Q(:,j)))
		enddo
		lam=1+lam/Tstar
		call lincg(Q,Qinv)
		if(IERCD()>0) call mdisp(real(lam).cvr.aimag(lam))
		a=1-[spread(lam,1,p1)]*[spread(lam,2,p1)]
		b=[spread(Qinv(:,p1),dim=2,ncopies=p1)*spread(Qinv(:,p1),dim=1,ncopies=p1)]
		val=reshape([matmul(Q,matmul(reshape(b/a,[p1,p1]),transpose(Q)))],[p1,p1])
	end function

	subroutine setstatespace1(th1,phi,theta,statevar)
		real			:: th1(nth1),Phi(p1,p1),theta(p1)
		real			:: poly(p1+1)
		real			:: statevar(p1,p1)
		integer		:: j
		
		Phi=0
		do j=2,p1
			Phi(j-1,j)=1
		enddo
		poly=realpolyfromroots(-getcfromh(th1(1:p1)))
		Phi(p1,:)=-poly(1:p1)
		if(p1>1) then
			theta=realpolyfromroots(-getcfromh(th1(p1+1:)))
		else
			theta=1
		endif
		Phi=eye(p1)+phi/Tstar
		statevar=getstatvar1(th1)
	end subroutine
	
	function getKalman_invar1(th1,y) result(val)
		real			:: th1(nth1),y(capN)
		real			:: val
		real			:: Phi(p1,p1),theta(p1)
		real			:: state(p1,2), statevar(p1,p1), vst(p1), cv, e(2)
		real			:: ssum(4)
		integer			:: s,j

		call setstatespace1(th1,phi,theta,statevar)
		state=0
		ssum=0
		j=1
		do s=1,Tstar
			state=matmul(Phi,state)
			statevar=matmul(Phi,matmul(statevar,transpose(Phi)))
			statevar(p1,p1)=statevar(p1,p1)+1.0
			if(s*capN>=j*Tstar) then
				vst=matmul(statevar,theta)
				cv=1.0/sum(vst*theta)
				e=[y(j),1.0]-matmul(theta,state)
				ssum(1:3)=ssum(1:3)+cv*[e(1)**2,e(2)**2,e(1)*e(2)]
				ssum(4)=ssum(4)+log(cv)
				state=state+matmul(reshape(vst,[p1,1]),reshape(cv*e,[1,2]))
				statevar=statevar-matmul(reshape(vst*cv,[p1,1]),reshape(vst,[1,p1]))
				j=j+1
			endif
		enddo
		val=.5*ssum(4)-0.5*log(ssum(2))-0.5*(capN-1)*log(ssum(1)-ssum(3)**2/ssum(2))
	end function
end module
	
	