#Assignment # 3
#For this assignment there are 3 different models:
#	-1) Schaefer production model with observation errors	(M1)
#	-2) Schaefer productio nmodel with process errors	(M2)
#	-3) Lagged Recruit Growth Survival Model			(M3)
#	-4) A total error version of the schaefer model 		(M4)

#rm(list=ls(all=T))		#will remove all objects from R-Workspace
setwd("~/Documents/UBC Courses/Fish 504/Tutorials2009")

#DATA SECTION
hake=read.table("Hakedata.txt",header=T)
ct = hake$Catch/1000
cpue=hake$CPUE
year=hake$Year
iyr=cbind(year,year+0.2,year+0.4,year+0.55,year+0.65)
L = 3

#PARAMETER SECTION
theta=c(log.k=log(2.9),log.r=log(0.38),log.q=log(0.35),tau=7)
phi=c(log.bo=log(3),log.k=log(11),logit.s=log(0.9/(1-0.9)),log.q=log(0.3),tau=10)
eta = list(log.k=log(2.9),log.r=log(0.38),log.q=log(0.35),tau=7,nu=rep(0,length=length(ct)))
i.eta=as.relistable(eta)

#PROCEDURE SECTION
M1=function(theta)
{
	with(as.list(theta),{
		bt=vector()
		k=exp(log.k); r=exp(log.r); q=exp(log.q)
		bt[1]=k; n=length(ct)
		for(i in 1:n)
		{
			bt[i+1]=bt[i]+r*bt[i]*(1-bt[i]/k)-ct[i]
		}
		epsilon=log(cpue)-log(q*bt[1:n])
		nloglike= -sum(dnorm(epsilon,0,1/tau,log=T))
		return(list(bt=bt,epsilon=epsilon,nloglike=nloglike,
			par=c(k=k,r=r,q=q,sig=1/tau)))
	})
}

M2=function(theta)
{
	with(as.list(theta),{
		bt=vector()
		k=exp(log.k); r=exp(log.r); q=exp(log.q)
		bt[1]=k; n=length(ct)
		for(i in 1:n)
		{
			bt[i+1]=cpue[i]/q+r*cpue[i]/q*(1-cpue[i]/q/k)-ct[i]
		}
		epsilon=log(cpue)-log(q*bt[1:n])
		nloglike= -sum(dnorm(epsilon,0,1/tau,log=T))
		return(list(bt=bt,epsilon=epsilon,nloglike=nloglike,
			par=c(k=k,r=r,q=q,sig=1/tau)))
	})
}


M3=function(theta)
{
	with(as.list(theta),{
		bt=rt=vector()
		bo=exp(log.bo); k=exp(log.k)+1.; s=1/(1+exp(-logit.s))
		a=k*(1-s); b=(k-1)/bo
		bt[1]=bo; n=length(ct)
		rt[1:L] = a*bo/(1+b*bo)
		for(i in 1:n)
		{
			if(i>L) rt[i]=a*bt[i-L]/(1+b*bt[i-L])
			bt[i+1]=s*bt[i]+rt[i]-ct[i]
		}
		if(exists("log.q")) {q=exp(log.q); epsilon=log(cpue)-log(q*bt[1:n])}
		if(!exists("log.q")){epsilon=log(cpue)-log(bt[1:n]); q=exp(mean(epsilon)); epsilon=epsilon-mean(epsilon)}
		if(exists("tau")) nloglike= -sum(dnorm(epsilon,0,1/tau,log=T))
		if(!exists("tau")) {tau = 1/sd(epsilon); nloglike= -sum(dnorm(epsilon,0,1/tau,log=T)) }
		return(list(bt=bt,epsilon=epsilon,nloglike=nloglike,
			par=c(bo=bo,k=k,s=s,q=q,sig=1/tau)))
	})
}

M4=function(theta)
{	#total errors (theta=(k,r,q,nu[1:n],tau,sig))
	with(relist(theta,skeleton=eta),{
		bt=vector()
		k=exp(log.k); r=exp(log.r); q=exp(log.q)
		bt[1]=k; n=length(ct)
		for(i in 1:n)
		{
			bt[i+1]=(bt[i]+r*bt[i]*(1-bt[i]/k)-ct[i])*exp(nu[i])
		}
		epsilon=log(cpue)-log(q*bt[1:n])
		nloglike= -sum(dnorm(epsilon,0,1/tau,log=T)) - sum(dnorm(nu,0,1/tau,log=T))
		return(list(bt=bt,epsilon=epsilon,nloglike=nloglike,
			par=list(k=k,r=r,q=q,sig=1/tau,nu=nu)))
	})
}


#function prototypes
fn1 = function(theta)		M1(theta)$nloglike
fn2 = function(theta)	M2(theta)$nloglike
fn3 = function(theta)	M3(theta)$nloglike
fn4 = function(theta)	M4(theta)$nloglike

#optimizations
if(!exists("f1")) f1 = optim(theta,fn1,method="BFGS",hessian=T)
if(!exists("f2")) f2 = optim(theta,fn2,method="BFGS",hessian=T)
if(!exists("f3")) f3 = optim(phi,fn3,method="BFGS",hessian=T)
if(!exists("f4")) f4 = optim(unlist(i.eta),fn4,method="BFGS",hessian=T)

#Question 1
print(round(M1(f1$par)$par,2));print(round(M2(f2$par)$par,2));print(round(M3(f3$par)$par,2))

#Question 2
epsilon=cbind(M1(f1$par)$epsilon,M2(f2$par)$epsilon,M3(f3$par)$epsilon)
matplot(iyr[,1:3],epsilon,type="h",xlab="Year",ylab="Residuals")
legend("topright",paste("std for M",1:3," = ",round(sd(epsilon),3),sep=""),bty="n",lty=1:3,col=1:3)

#Question 3
aic = c(2*M1(f1$par)$nloglike+2*length(f1$par),
	    2*M2(f2$par)$nloglike+2*length(f1$par),
	    2*M3(f3$par)$nloglike+2*length(f3$par)+1)  #+1 because L is also a parameter
legend("bottomright",paste("AIC for M",1:3," = ",round(aic,3),sep=""),bty="n",lty=1:3,col=1:3)

#Question 4 (a few ways to do this but M4 uses a total error approach where variance in 
#observation errors is assumed to equal variance in process errors).
epsilon=cbind(M1(f1$par)$epsilon,M2(f2$par)$epsilon,M3(f3$par)$epsilon,M4(f4$par)$epsilon,M4(f4$par)$par$nu)
matplot(iyr,epsilon,type="h",xlab="Year",ylab="Observation error residuals")
legend("topright",paste("std for M",c(1:4,4)," = ",round(sd(epsilon),3),sep=""),bty="n",lty=1:5,col=1:5)


#So you want to do likelihood profiles for MSY in the LRGS model...Here's how I did it.
#Note that this method is very different from what is used in the Ecological Detective (ED).
#I've parameterized the model interms of MSY and Fmsy, then derive Bo and kappa given s.
#For each MSY hypothesis I then calculate the mle estimates of Fmsy and s using optim and 
#use the conditional MLE for q and Tau (nuisance parameters).
#Whereas, in the ED, the estimate Bo,k,s using optim ad derive MSY given Bo,k,s and add a
#penalty (p) to the negative loglikelihood p=k(MSY(Bo,k,s)-MSY)^2, where k is a large penalty
#weight to ensure find Bo,k,s parameters that explain the data and the MSY hypothesis.

get.k.bo<-function(theta)
{	#This function returns the bo, and k values given
	#initial estimates of MSY and Fmsy, and s
	#Parameter vector is: c(MSY,Fmsy,logit.s,log.q,tau)
	with(as.list(theta),{
		s=1/(1+exp(-logit.s))
		a = -(1-s+Fmsy)^2/(-1+s)
		b = -Fmsy*(1-s-a+Fmsy)/(MSY*(1-s+Fmsy))
		bo= -(-1+s+a)/(b*(-1+s))
		k = a/(1-s)
		return(c(bo,k))
	})

}


iphi=c(MSY=0.26,Fmsy=0.3,logit.s=log(0.9/(1-0.9)),log.q=log(0.3),tau=10)
mop=function(theta){
	#This function first gets bo and k, given MSY and FMY and s (calls get.k.bo(theta))
	#then constructs a parameter vector phi, given the bo and k values
	#next, call the M3 model with phi as the argument.
	with(as.list(theta),{
		bo.k=get.k.bo(theta)
		phi=c(log.bo=log(bo.k[1]),log.k=log(bo.k[2]),logit.s=logit.s,log.q=log.q,tau=tau)
		#add a prior here for Fmsy.
		fprior= -dbeta(Fmsy,1.001,1.001,log=T)  #without this prior/penalty Fmsy>1
		M=M3(phi); M$nloglike=M$nloglike+fprior
		return(M)
	})
}
fnx=function(theta) mop(theta)$nloglike  #prototype for the negative loglike
f5=optim(iphi,fnx,method="BFGS",control=list(maxit=500),hessian=T)
V=solve(f5$hessian)
std=sqrt(diag(V))
R=V/(std%o%std)
print(round(R,4))

msy = seq(0.26,0.4,length=50)
msy.profile=function(msy)
{
	fn=function(theta){
		with(as.list(theta),{
			Fmsy=1/(1+exp(-logit.Fmsy))
			iphi=c(MSY=msy,Fmsy=Fmsy,logit.s=logit.s)
			bo.k=get.k.bo(iphi)
			phi=c(log.bo=log(bo.k[1]),log.k=log(bo.k[2]),logit.s=logit.s)
			#add a prior here for Fmsy.
			fprior= -dbeta(Fmsy,1.001,1.001,log=T)  #without this prior/penalty Fmsy>1
			#Ensure MSY likelihood profile is not sensitive to the fprior.
			M=M3(phi); M$nloglike=M$nloglike+fprior
			return(M$nloglike)
		})
	}
	iphi=c(logit.Fmsy=log(0.95/(1-0.95)),logit.s=log(0.86/(1-0.86)))#,log.q=log(0.3),tau=10)
	fit=optim(iphi,fn,method="BFGS",control=list(maxit=500),hessian=F)
	print(c(msy,fit$convergence,fit$counts))
	return(fit$value)
}
lvec=sapply(msy,msy.profile)
plot(msy,exp(-(lvec-min(lvec))),type="l")
plot(msy,pchisq(2*(lvec-min(lvec)),df=1),type="l",main="Chi-square probability")