module basicmodule
	use myfuncs
	use dotops
	use ML

	logical, parameter	:: CYrun=.false.		! flag for Campbell Yogo application
	integer, parameter	:: Tstar=1000
	integer				:: capN
	integer				:: p
	integer, parameter	:: pmax=5
	real				:: hlb=0.0, hub		! lower and upper bound on each element in h
	real, parameter		:: hhlb=3, hhub=50	! lower and upper bound on half-life
	real				:: hhf=1			! for quarter-life, set hhf=2, for sqrt(0.5)-life, hhf=0.5
	integer				:: hhspecialflag=0	! 0 default, 1 use sum of AR coefficients, 2 absolute sum of correlations
	real				:: maxhhfac			! how many multiples of Tstar do we compute correlations to obtain half-lives
	real				:: priorkappa		! smoothness parameter for baseline prior
	logical				:: rawpriorflag		! if true, generate draws from baseline prior (only used to generate Figure 3)
	
	integer				:: icase

	real, allocatable	:: ydata(:)
	integer				:: Tdata
	
	integer, parameter	:: nhhbins=39		! bins for adjusted prior so that implied prior on half-life is uniform
	real, parameter		:: hhgrid(nhhbins+1)=[(hhlb+i*(hhub-hhlb)/nhhbins,i=0,nhhbins)]
	real				:: hhprior0(nhhbins,pmax)
	logical				:: hhpriorset
	
	integer, parameter	:: ns=10000, nruns=20,nskip=10	! number of posterior draws and number of chains
	integer, parameter	:: nsa=ns*nruns
	
	character(len=100)	:: dir="c:/dropbox/mystuff/beyondltu/replication/"
	
	integer				:: acc
	real				:: stepsize
!$omp threadprivate(acc,stepsize,p)
	
	type tpara 
		real	:: hAR(pmax),hMA(pmax-1)		! h-parameters for GLTU model
	end type
	
	contains

	function getcfromh(h) result(c)		! get c (or g) from h
		implicit none
		real		:: h(:)
		complex	:: c(size(h))
		real		:: h1,h2
		complex	:: d
		integer	:: i
		do i=1,size(h)/2
			h1=h(2*i-1)
			h2=h(2*i)
			d=sqrt(cmplx(h2**2 - h1**2))
			c(2*i-1:2*i)=[h2-d,h2+d]
		enddo
		if(mod(size(h),2)==1) then
			c(size(h))=h(size(h))
		endif
	end function
	
end module
	

module iomodule
!DIR$ NOOPTIMIZE
	use basicmodule
	implicit none
	
	contains
	
	subroutine loaddata
		real, allocatable	:: mdata(:)
		integer	:: s,j,T
		if(CYrun) then
			mdata=loadvec(trim(dir)//"data/price_dividend.txt")
		else		
			mdata=loadvec(trim(dir)//"data/USUK_RER_long.txt")
		endif
		T=size(mdata)
		Tdata=T
		if(allocated(ydata)) deallocate(ydata)
		allocate(ydata(capN))
		print *,"length of data in years:", Tdata
		j=1
		do s=1,T
			if(s*capN>=j*T) then
				ydata(j)=mdata(s)
				j=j+1
			endif
		enddo
		call mdisp(mdata)
		call mdisp(ydata)
	end subroutine
	
	function kerneldens(x,draws,bw,xmin,xmax) result(val)
		real		:: x,draws(:),bw,xmin,xmax,val
		val=sum(gausspdf((x-draws)/bw))/(size(draws)*bw)
		val=val/(gausscdf((xmax-x)/bw)-gausscdf((xmin-x)/bw))
	end function
	
	subroutine mkpostplot(hhdraws)
		real		:: hhdraws(:)
		integer	:: i
		integer, parameter	:: npoints=200
		real	,save	:: plot(0:5,0:npoints)

		do i=0,npoints
			plot(0,i)=hhlb+i*(hhub-hhlb)/npoints
			plot(p,i)=kerneldens(plot(0,i),hhdraws,1.0,hhlb,hhub)
		enddo
		print *,"posterior plot for p=",p
		call mdisp(plot(0:p,:))
		call storeML(plot(0,:),"x")
		call storeML(plot(1:p,:),"y")
		call execinML("figure(11)")
		call execinML("plot(x,y)")
		if(p==pmax) then 
			call savemat(trim(dir)//"postplot"//convtos(icase)//".txt",plot)
		endif
	end subroutine
	
	subroutine plotspecdens(paras)		! visualize implied spectral densities
		type(tpara)	:: paras(:)
		integer	:: id
		integer, parameter	:: npoints=200
		real, parameter		:: lammax=200
		integer				:: i
		real, parameter		:: lamgrid(npoints)=[(lammax*(i-1.0)/npoints,i=1,npoints)]
		real		:: mp(npoints,size(paras))
		integer	:: n
		
		n=size(paras)
		do i=1,n
			mp(:,i)=getldens(paras(i))
		enddo
		call storeML(mp,"y")
		call storeML(lamgrid,"x")
		call execinML("plot(x,y)")
	contains
	
		function getldens(para) result(val)
			use polymod
			type(tpara)	:: para
			real		:: val(npoints)
			real		:: apoly(p+1), bpoly(p), aapoly(2*p+1), bbpoly(2*p-1)
			integer		:: i
			apoly=realpolyfromroots(-getcfromh(para%hAR(1:p)))
			aapoly=polymult(apoly,polyminus(apoly))
			if(p==1) then
				bpoly=1
				bbpoly=1
			else
				bpoly=realpolyfromroots(-getcfromh(para%hMA(1:p-1)))
				bbpoly=polymult(bpoly,polyminus(bpoly))
			endif
			do i=1,npoints
				val(i)=polyevalcmplx(bbpoly,cmplx(0,lamgrid(i)))/polyevalcmplx(aapoly,cmplx(0,lamgrid(i)))
			enddo
			val=log(val)
		end function
	end subroutine
			
end module
