#This routine calculates some basic risk metrics using the DennisHolmes method (aka Holmes and Fagan 2002)
#Confidence intervals on mu, Probability of extinction with 95% confidence intervals,
#Probability of decline with 95% confidence intervals, Probability of observing extinction, 
#Probability of observing a given decline
#
#mu: estimate of parameter mu
#sigma2: estimate of parameter sigma2_p=sigma2_slp
#sigma2np: estimate of paramter sigma2_np
#nq: current population size; absolute for extinction calcs or relative for decline calcs
#ne: forecast population size; absolute for extinction calcs or relative for decline calcs
#n: number of counts
#L: filter length
#years: number of forecast years; used in probability of decline and extinction calcs
#probofext: F means do not do prob of extinction calculations
#probofdecline: F means do not do prob of decline calculations
#
#Example 1: To calculate the probability of 90% decline in 25 years for a population 
#   riskcalcs(mu=-.05, sigma2=.05, nq=10, ne=1, n=20, L=4, years=25, probofdecline=T)
#Example 2: To calculate the probability of extinction in 25 years for a population of size 200
#   riskcalcs(mu=-.05, sigma2=.05, nq=200, ne=1, n=20, L=4, years=25, probofext=T)
#Example 3: To calculate the probability of OBSERVING extinction in 25 years for a population of size 200
#   riskcalcs(mu=-.05, sigma2=.05, sigma2np=.2, nq=200, ne=1, n=20, L=4, years=25, probofext=T)

riskcalcs_function(mu=NA, sigma2=NA, sigma2np=0, nq=NA, ne=NA, 
	n=NA, L=NA, years=1, probofext=F, probofdecline=F)
{
	if(is.na(mu)){print("Error denniscalcs; need mu"); break}
	if(is.na(sigma2)){print("Error denniscalcs; need sigma2"); break}
	if(is.na(nq)){print("Error denniscalcs; need nq"); break}
	if(is.na(ne)){print("Error denniscalcs; need ne"); break}
	if(is.na(n) | is.na(L)){print("Error denniscalcs; need n and L"); break}
	if(n<10) {print("Error denniscalcs; n must be greater than 10"); break}
		
	if(sigma2==0) sigma2_.000001 #will choke on 0
	tq_n-L+1
	if(n==10) dfRslp_1.42 
		else dfslp_.333+.212*n-.387*L  #from Holmes and Fagan 2002
	q _ tq-1;
	d _ .000001	
	cumprobext_NA
	cumprobextwithobserror_NA
	cumprobdec_NA
	cumprobdecwithobserror_NA
	cumprobext95confints_NA
	cumprobdec95confints_NA
	mu95confints_NA
	
	#Get confidence intervals on lambda
		muup95 _ ( mu - qt(.025,dfslp)*sqrt(sigma2*(1/q)) )
		mulow95 _ ( mu + qt(.025,dfslp)*sqrt(sigma2*(1/q)) )
		mu95confints _ c(muup95,mulow95)
	
	#Get extinction probabilities on true population size
	if(probofext)
	{
		tmp_2*log(nq/ne)*abs(mu)/sigma2
		tmp_tmp/700 #this is a trick to get rid of NAs when tmp is big
		if(tmp<1)tmp_1
		G _ pnorm( (-log(nq/ne) + abs(mu)*years)/sqrt(sigma2*years)) + (exp((1/tmp)*2*log(nq/ne)*abs(mu)/sigma2)*
		pnorm( (-log(nq/ne) - abs(mu)*years)/sqrt(sigma2*years))^(1/tmp))^tmp;
		if(mu < 0) probext _1 else probext _ (ne/nq)^(2*mu/(sigma2*(q-1)/q))
		cumprobext _ G*probext;
	}
	
	#Get probability true population is nq/ne lower at time years
	if (probofdecline) #Prob that stock is nq/ne lower at year t
		{
		if(sigma2np < 0) sigma2np_0
		G _ 1-pnorm( (log(nq/ne) + mu*years)/sqrt(sigma2*years));
		cumprobdec _ G;
		}

	#Get confidence intervals on prob of decline or extinction via monte carlo simulation following dennis 1991
	if (probofext | probofdecline)
	{
		tmpreps_1000
		murnds_rt(tmpreps,dfslp)*sqrt(sigma2/(n-L))+mu
		muhigh_qt(.975,dfslp)*sqrt(sigma2/(n-L))+mu
		mulow_qt(.025,dfslp)*sqrt(sigma2/(n-L))+mu
		s2slprnds_sigma2*dfslp/rchisq(tmpreps,dfslp)
		s2low_sigma2*dfslp/qchisq(.975,dfslp)
		s2high_sigma2*dfslp/qchisq(.025,dfslp)
		okparams_(murnds<=muhigh) & (murnds>=mulow) & (s2slprnds<=s2high) & (s2slprnds>=s2low)
		murnds_murnds[okparams]
		s2slprnds_s2slprnds[okparams]
		if(probofext)
			{
			cumprobext95confints_{}
			for (i in 1:length(years))
				{
				tmp_2*log(nq/ne)*abs(murnds)/s2slprnds
				tmp_tmp/700 #this is a trick to get rid of NAs when tmp is big
				tmp[tmp<1]_1
				G _ pnorm( (-log(nq/ne) + abs(murnds)*years[i])/sqrt(s2slprnds*years[i])) + (exp((1/tmp)*2*log(nq/ne)*abs(murnds)/s2slprnds)*
					pnorm( (-log(nq/ne) - abs(murnds)*years[i])/sqrt(s2slprnds*years[i]))^(1/tmp))^tmp;
				probext _ (ne/nq)^(2*murnds/(s2slprnds*(q-1)/q))
				probext[murnds<0]_1
				cumprobextrnds _ G*probext;
				upper95 _ quantile(cumprobextrnds,.975);
				low95 _ quantile(cumprobextrnds,.025);
				cumprobext95confints_c(cumprobext95confints,quantile(cumprobextrnds,c(.975,.025)))
				}
			}#if probofext=T
		if(probofdecline)
			{
			cumprobdec95confints_{}
			for (i in 1:length(years))
				{
				cumprobdecrnds _ 1-pnorm( (log(nq/ne) + murnds*years[i])/sqrt(s2slprnds*years[i]));
				cumprobdec95confints_c(cumprobdec95confints,quantile(cumprobdecrnds,c(.975,.025)))
				}
			}#if probofdecline=T

	}	#probconfints
	

	#Probability of observing the stock is below threshold with time tau given that observations have error
	if (probofext==T & sigma2np!=0) 
	{	
		prob_0
		reps_500
		tau_1/1 #basically the step length
		projtime_years/tau
		for(nn in 1:reps)
			{
			N_rep(0,projtime+1)
			N[1]_100
			ep_exp(rnorm(projtime+1,mu*tau,sqrt(sigma2*tau)))
			if(sigma2np<=0) enp_1 else enp_exp(rnorm(projtime+1,0,sqrt(sigma2np)))
			for(year in 2:(projtime+1))
				N[year]_N[year-1]*ep[year]
			N_enp*N
			if(min(N[seq(1/tau,projtime,1/tau)])<=(ne/nq)*N[1]) prob_prob+1
			}
		cumprobextwithobserror_prob/reps
	}
	
	#Prob of observing stock below threshold at time tau given that observations are corrupted
	if (probofdecline==T & sigma2np!=0) #Prob that stock is nq/ne lower at year t
	{	
		reps_500
		projtime_years
		prob_0
		for(nn in 1:reps)
			{
			N_rep(0,projtime+1)
			N[1]_100
			ep_exp(rnorm(projtime+1,mu,sqrt(sigma2)))
			if(sigma2np<=0) enp_1 else enp_exp(rnorm(projtime+1,0,sqrt(sigma2np)))
			for(year in 2:(projtime+1))
				N[year]_N[year-1]*ep[year]
			N_enp*N
			if(N[projtime+1]<=(ne/nq)*N[1]) prob_prob+1
			}
		cumprobdecwithobserror_prob/reps
	}

	output _ list(mu95confints=mu95confints, cumprobext=cumprobext, cumprobext95confints=cumprobext95confints, 
		cumprobdec=cumprobdec, cumprobdec95confints=cumprobdec95confints,
		cumprobextwithobserror=cumprobextwithobserror,
		cumprobdecwithobserror=cumprobdecwithobserror	)
	
	return(output)
}