# This function estimates parameters assuming no inputs into the population
# No missing years allowed in N; Interpolate if there are missing years

DHMethod_function(N, filtlen=4) #default filter length is 4
{
####Check if there is enough data; break if not
if(length(N)<(filtlen+8))
	{
	print("Error: not enough data");
	break
	}

####Set up running sum filter and runsum
runsumfilt _ rep(1,filtlen)
runsum _ filter(N,runsumfilt,method="con",sides=1);
runsum_runsum[is.number(runsum)] #filter func puts in trailing NAs; remove
runsum_as.vector(runsum) #remove Splus timeseries data format
if(is.element(0,runsum))
	{print("Error: Zeros in runsum"); break} 

####Get logratios of N and runsum at tau=1
runlen_length(runsum)
Rlogratios _ log(runsum[2:runlen])-log(runsum[1:(runlen-1)]);
lenN_length(N)
Nlogratios_log(N[2:lenN]/N[1:(lenN-1)])
Nlogratiosyrs_(1:length(Nlogratios))[is.finite(Nlogratios)] #indices where Nlogratios is not zero; need for trend test
Nlogratios_Nlogratios[is.finite(Nlogratios)] #chuck cases where N is zero

	
####Do diagnostics on data
#check for normality
	normfit_ks.gof(Nlogratios,dist="norm")
	print(paste("p.value for normality goodness of fit test:", 
		format(normfit$p.value,digits=4)))

#Test for trends in data
	trendstat_ls.print(lsfit(Nlogratiosyrs,Nlogratios),print.it=F)
	trendpvalue_trendstat$summary[7]
	names(trendpvalue)_NULL
	print(paste("p.value for trend test:",format(trendpvalue,digits=4)))

#check for 1st order autocorrelation; detrend data before doing this
	lenNlogs_length(Nlogratios)
	logratiodetrend_
		Nlogratios-(trendstat$coef.table[1,1]+trendstat$coef.table[2,1]*
		(1:lenNlogs))
	tmpacf_cor.test(logratiodetrend[1:(lenNlogs-1)], 
		logratiodetrend[2:lenNlogs],method="spearman")
	print(paste("p.value for autocorrelation test:",
		format(tmpacf$p.value,digits=4)))

#check for and remove outliers in runsum logratios which is used to calc muR
	mulsfit_lsfit(rep(1,(runlen-1)),Rlogratios,intercept=F)
	mulsdiag_ls.diag(mulsfit)
	outliercutoff_qt(.9,runlen-1) #this is a generous; qt(.975) is standard
	mudfits_abs(mulsdiag$dfits)<outliercutoff 
	Rlogratios_Rlogratios[mudfits]
	if(prod(mudfits)==0) print("Outliers in running sum logratios")

####Estimate mu, sigma2p, sigma2np, and sigma2muR
#get mu
	muR_mean(Rlogratios)

#Get sigma2p using slope of variance versus tau	
	taumax_4
	runlen_length(runsum)
	varrunsums_c()
	for (i in 1:taumax)
		{
		tmp _ log(runsum[(i+1):runlen])-log(runsum[1:(runlen-i)]);
		varrunsums _ c(varrunsums,var(tmp));
		}
	#Get variance slope by linear regression
	lmsigma2 _ lm(varrunsums~c(1:length(varrunsums)))
	sigma2slp_lmsigma2$coefficients[2];
	names(sigma2slp)_NULL #cleans up presentation

	#comment if the sigma regression is bad
	lmsigma2sum_summary(lmsigma2);
	print(paste("r-squared for linearity of variance slope:",
		format(lmsigma2sum$r.squared,digits=4)))

	#If sigma2slp is basically 0, can get sigma2<0 if answer should really
	#be sigma2=0; check and replace if applicable
	comment_""
	if (sigma2slp < 0 & lmsigma2sum$sigma > .01) 
		print("Sigma2slp less than 0 and not zero"); 
	if (sigma2slp < 0 & lmsigma2sum$sigma < .01) 
		{
		sigma2slp_0.00001; 
		print("Replaced neg sigma2slp with 1E-5")
		}

#Get sigma2np
	sigma2np_(var(Nlogratios)-sigma2slp)/2

#Get sigma2muR
	sigma2muR_(2*sigma2np/filtlen + (lenN-filtlen)*sigma2slp)/(lenN-filtlen)^2
	
return(list(muR=muR, sigma2slp=sigma2slp, sigma2np=sigma2np, 
	sigma2muR=sigma2muR, n=lenN, L=filtlen))
}