*! spectdens computes periodogram and spectral density estimate
*! v.1.0 Jan 8, 2016 Huseyin Tastan (tastan@yildiz.edu.tr)
*! v.1.1 May 25, 2016 Fixed the problem with saving output dta file

cap program drop spectdens
program define spectdens, rclass
	version 12
	syntax varname(ts) [if] [in] [, 								///
									Meanadj 						///
									DETrend 						///
									Weights(string) 				///
									Kernel(string) 					///
									Bandwidth(integer 0) 			///
									CONVdaniell(numlist >0 integer) ///
									CInterval 						///
									log 							///
									out(string) 					///
									noGRAPH                         ///
									replace ]
	marksample touse
	qui count if `touse'
	local n = r(N)
	local m = floor(`n'/2)
	local mm = `m' + 1
	if `mm'>400 {
		set matsize `mm'
		}
	if `n' == 0 {
		error 2000
	}
	_ts timevar, sort
	markout `touse' `timevar'
	tsreport if `touse', report
	if r(N_gaps) {
		di in red "Sample may not contain gaps"
		exit
	}


local nogr "`graph'"  // default is graph
local ww "`weights'"

if "`ww'" != "" & "`kernel'" != "" & "`convdaniell'" != ""{
	di in red "Choose only one: weights or kernel type-bandwidth or convolution"
	error 184
}
if "`ww'" != "" & "`kernel'" != "" {
	di in red "Choose only one: weights or kernel type-bandwidth or convolution"
	error 184
}
if "`kernel'" != "" & "`convdaniell'" != ""{
	di in red "kernel type-bandwidth cannot be used with convolution"
	error 184
}
if "`ww'" != "" & "`convdaniell'" != ""{
	di in red "Weights and convdaniell cannot be used together"
	error 184
}
if "`ww'" == "" & "`kernel'" == "" & "`convdaniell'" == "" {
	local kernel "daniell" // default is Daniell kernel if no option input
	di "Using default Daniell (uniform) kernel"
}
if "`ww'" != "" & "`kernel'" == "" & "`convdaniell'" == ""{
* check if weights supplied by the user include < or >
	if regexm("`ww'","<")==0 {
		di as err "missing < in weights()"
		di as err "weights(`weights') invalid"
	}
	if regexm("`ww'",">")==0 {
		di as err "missing > in weights()"
		di as err "weights(`weights') invalid"
	}
	if regexm("`ww'","/")==1 {
		di as err "Do not use / in weights()"
		di as err "weights(`weights') invalid"
		error 184
	}
	local nw = wordcount("`ww'")
	if mod(`nw',2)!=1 {
		di as err "Total window length in weights() should be an odd number"
		di as err "weights(`weights') invalid"
	}

	local wei = regexr(regexr("`ww'", "<", " "), ">", " ")
	tempname wvec
	mat `wvec' = J(`nw',1,0)
	local bw = (`nw'-1)/2
	local ind=0
	foreach v of local wei {
		local ind = `ind'+1
		mat `wvec'[`ind',1] = `v'
	}
local nte "Weigths = `ww'"
}


if "`convdaniell'" != "" & "`ww'" == "" & "`kernel'" == "" {
	tempname conv wvec
	local np: word count `convdaniell'
	tokenize `convdaniell'
		if `np'<1 {
		di in red "Convolution incorrectly specified"
		error 184
		}
		if `np'==1 {
		di "Using Daniell kernel with window `1'"
		_convdaniell `1' 0
		mat `conv' = c
		}
		if `np'==2 {
		_convdaniell `1' `2'
		mat `conv' = c
		}
		if `np'>2 {
			forv k=1(1)`np' {
				tempname u`k' v`k'
				_convdaniell ``k'' 0
				mat `u`k'' = c
			}
			_convvector `u1' `u2'
			mat `v1' = c
			forv k=3(1)`np' {
				local kk = `k'-1
				local kkk = `k'-2
				_convvector `v`kkk'' `u`k''
				mat `v`kk'' = c
			}
		mat `conv' = `v`kk''
		}

mat `wvec' = `conv'
local nte "Convolution of Daniell kernel, lags = `convdaniell'"
return local kernel `"convdaniell"'
return local lags = `"`convdaniell'"'

}

if "`kernel'"!="" & "`ww'" == "" & "`convdaniell'" == ""{
	if "`kernel'"=="daniell" {
		local bw = ceil(0.25*`mm'^(1/5))
		if `bandwidth' !=0 {
			local bw = `bandwidth'
		}

		tempname wvec
		mat `wvec' = J((2*`bw'+1),1,0)
		local ind=0
		forv q =-`bw'(1)`bw' {
				local bb = 1/(2*`bw'+1)
				local ind = `ind'+1
				mat `wvec'[`ind',1] = `bb'
		}
	}
	if "`kernel'"=="mdaniell" {
		local bw = ceil(0.25*`mm'^(1/5))
		if `bandwidth' !=0 {
			local bw = `bandwidth'
		}
		tempname wvec
		mat `wvec' = J((2*`bw'+1),1,0)
		local ind=0
		forv q =-`bw'(1)`bw' {
			local ind = `ind'+1
			if `q'==-`bw'{
				mat `wvec'[`ind',1] = 1/(4*`bw')
			}
			if `q'>-`bw' & `q'<0{
				mat `wvec'[`ind',1] = 1/(2*`bw')
			}
			if `q'==0{
				mat `wvec'[`ind',1] = 1/(2*`bw')
			}
			if `q'<`bw' & `q'>0{
				mat `wvec'[`ind',1] = 1/(2*`bw')
			}
			if `q'==`bw'{
				mat `wvec'[`ind',1] = 1/(4*`bw')
			}
		}
	}
	if "`kernel'"=="bartlett" {
		local bw = ceil(0.5*`mm'^(1/3))
		if `bandwidth' !=0 {
			local bw = `bandwidth'
		}
		tempname wvec
		mat `wvec' = J((2*`bw'+1),1,0)
		local ind=0
		forv q =-`bw'(1)`bw' {
			local ind = `ind'+1
			mat `wvec'[`ind',1] = 1-abs(`q'/(`bw'+1))
		}
	}
	if "`kernel'"=="parzen" {
		local bw = ceil(`mm'^(1/5))
		if `bandwidth' !=0 {
			local bw = `bandwidth'
		}
		tempname wvec
		mat `wvec' = J((2*`bw'+1),1,0)
		local ind=0
		forv q =-`bw'(1)`bw' {
			local ind = `ind'+1
			local c = abs(`q'/(`bw'+1))
				if `c'>=0 & `c'<0.5 {
					local bb = 1-6*`c'^2+6*`c'^3
				}
				else if `c'>=0.5 & `c'<=1 {
					local bb = 2*(1-`c')^3
				}
				else {
					local bb = 0
				}
			mat `wvec'[`ind',1] =  `bb'
		}
	}
	if "`kernel'"=="tukey" {
		local bw = ceil((2/3)*`mm'^(1/5))
		if `bandwidth' !=0 {
			local bw = `bandwidth'
		}
		tempname wvec
		mat `wvec' = J((2*`bw'+1),1,0)
		local ind=0
		forv q =-`bw'(1)`bw' {
			local ind = `ind'+1
			local c = `q'/(`bw'+1)
				if `c'<=1 {
					local bb = (1+cos(_pi*`c'))/2
				}
				else {
					local bb = 0
				}

			mat `wvec'[`ind',1] =  `bb'
		}
	}
	if "`kernel'"=="qs" {
		local bw = ceil((1/2)*`mm'^(1/5))
		if `bandwidth' !=0 {
			local bw = `bandwidth'
		}
		tempname wvec
		mat `wvec' = J((2*`bw'+1),1,0)
		local ind=0
		forv q =-`bw'(1)`bw' {
			local ind = `ind'+1
			local c = abs(`q'/`bw')
				if `c'!=0 {
					local bb = (25/(12*_pi^2*`c'^2))*((sin(6*_pi*`c'/5)/(6*_pi*`c'/5)) - cos(6*_pi*`c'/5))
				}
				else {
					local bb = 1
				}
			mat `wvec'[`ind',1] =  `bb'
		}
	}
local nte "kernel(lag) = `kernel'(`bw')"
return local kernel `"`kernel'"'
return scalar halfbw = `bw'
}

// change row names of the weight vector
local rnames = ""
local rr=(rowsof(`wvec')-1)/2
forv q =-`rr'(1)`rr' {
	local rnames = "`rnames' `q'"
}
mat rownames `wvec' = `rnames'

preserve

tempvar x xvar trend freq omega wave k xa xb Preal Pimaginary P Pgraph P1 S ///
		S1 SL SU logS logP CI
tempname Pmat Smat

qui{
	keep if `touse'
	gen double `x' = `varlist'

	if "`meanadj'" == "" & "`detrend'" == "" {
		gen double `xvar' = `x' if `touse'
	}
	if "`meanadj'" != "" & "`detrend'" == "" {
		sum `x' if `touse', meanonly
		gen double `xvar' = `x'-r(mean) if `touse'
	}
	if "`detrend'" != "" {
		gen double `trend' = _n if `touse'
		reg `x' `trend' if `touse'
		predict `xvar' if `touse', resid
	}
	gen double `k' = _n - 1 in 1/`mm'
	gen double `freq' = `k'/`n'
	gen double `omega' = 2*_pi*`freq'
	gen double `wave' = 1/`freq'
	label var `omega' "Frequency"
	label var `wave' "Wavelength (Period)"
	label var `k' "k=0,1,2,...,floor(N/2)"
    label var `omega' "Fourier Frequencies"
	label var `freq' "Natural Frequencies"

	fft `xvar' if `touse', gen(`xa' `xb')
	gen double `Preal' = (2/`n')*`xa' in 1/`mm'
	gen double `Pimaginary' = -(2/`n')*`xb' in 1/`mm'
	label var `Preal' "Real part of periodogram"
	label var `Pimaginary' "Imaginary part of peridogram"

	gen double `P' = (`n'/2)*(`Preal'^2 + `Pimaginary'^2)
	local P1 = `P'[1]
	replace `P' = `P'[2] in 1 /* use second pergram ordinate in place of first*/
	                          /* (omega=0) for smoothing */
	tsset `k'
	keep in 1/`mm'
	gen double `S' = .
	* use cyclical endpoints
	mata: smoothdens("`P'", "`S'")
	local wssq = r(wssq)
	local dof = 2/`wssq'

	replace `P' = `P1' in 1 // re-replace periodogram at freq 0
	gen double `Pgraph' = `P'/(4*_pi) // scale periodogram ordinate for graphing
	replace `Pgraph' = . in 1 // re-replace a missing value at freq 0 for graphing
	label var `P' "Periodogram"
	label var `S' "Spectral Density"
	label var `Pgraph' "Periodogram"

	* in logs
	if "`log'" == ""{
		gen double `SL' = `dof'*`S'/invchi2(`dof',0.975)
		gen double `SU' = `dof'*`S'/invchi2(`dof',0.025)
	}
	else {
		gen double `logS' = log(`S')
		gen double `SL' = `logS'+ log(`dof'/invchi2(`dof',0.975))
		gen double `SU' = `logS'+ log(`dof'/invchi2(`dof',0.025))
		gen double `logP' = log(`P')
	}
}

if "`nogr'" == ""{
	if "`cinterval'" == "" & "`log'" == ""{
		graph twoway  (scatter `Pgraph' `freq',         		  ///
									    msymbol(smcircle)         ///
										mcolor(maroon)            ///
										xaxis(1 2) yaxis(1 2))    ///
					   (line `S' `freq', 						  ///
										lwidth(medthick) 		  ///
										lcolor(blue)   			  ///
										xaxis(1 2) yaxis(1 2))    ///
                   , xmtick(#25) scale(0.9)                       ///
				     xtitle("Natural Frequencies", axis(1))       ///
				     xtitle("", axis(2)) legend(order(2 1))		  ///
					 ytitle("Spectral Density", axis(1))          ///
					 graphregion(color(white)) 					  ///
					 title("Spectral Density of `varlist'")       ///
					 note("`nte'")
	}
	else if "`cinterval'" == "" & "`log'" != ""{
		graph twoway (line `logS' `freq', xaxis(1 2)	 		     ///
							lcolor(blue) 					 		 ///
							yaxis(1 2)						 		 ///
							lwidth(medthick) 				 		 ///
							xtick(#25) 						 		 ///
							ytitle("Log Spectral Density")   		 ///
							ytitle("", axis(2))              		 ///
						, scale(0.9) xtitle("", axis(2)) 		     ///
						  legend(cols(3))                  		     ///
						  graphregion(color(white)) 		 	     ///
						  xmtick(#25)                      		     ///
						  title("Log-Spectral Density of `varlist'") ///
						  note("`nte'"))
	}
	else if "`cinterval'" != "" & "`log'" != ""{
		graph twoway (rarea `SL' `SU' `freq', 							 ///
							astyle(ci) 							    	 ///
							lcolor(black) 						    	 ///
							legend(label(1 "95% CI"))       	    	 ///
							legend(label(2 "Log Spectral Density"))		 ///
							xaxis(1 2) yaxis(1 2))    				  	 ///
					 (line `logS' `freq', 								 ///
							lcolor(blue) 								 ///
							lwidth(medthick) 						     ///
							xaxis(1 2) yaxis(1 2)    					 ///
							ytitle("Log of Smoothed Spectral Density"))  ///
				     , legend(order(2 1)) legend(cols(3)) 				 ///
						xtick(#25) scale(0.9)							 ///
						note("`nte'")		 							 ///
						graphregion(color(white)) 						 ///
						xtitle("", axis(2)) 							 ///
						title("Log-Spectral Density of `varlist'")
	}
	else {
		graph twoway (rarea `SL' `SU' `freq', xaxis(1 2) yaxis(1 2)	   	 ///
											  astyle(ci) 		 		 ///
											  lcolor(black) 	 		 ///
											  legend(label(1 "95% CI"))) ///
					 (scatter `Pgraph' `freq',         		  		     ///
											   msymbol(smcircle)  		 ///
											   mcolor(maroon)    		 ///
											   xaxis(1 2) yaxis(1 2))    ///
					 (line  `S' `freq', 						  		 ///
										lwidth(medthick) 		  		 ///
										lcolor(blue)   			  		 ///
										xaxis(1 2) yaxis(1 2))        	 ///
                     , scale(0.9) xtitle("Natural Frequencies", axis(1)) ///
					   xtitle("", axis(2)) 								 ///
					   ytitle("Spectral Density", axis(1)) 			     ///
					   xtick(#25)										 ///
				       graphregion(color(white)) 						 ///
					   legend(order(3 2 1)) legend(cols(3)) 			 ///
					   note("`nte'")		 							 ///
					   title("Spectral Density of `varlist'")
	}
}
else if "`nogr'" != ""{
	di "No graphical output"
	di as text "{hline 72}"
}



if "`log'" == ""{
	mkmat `freq' `S', matrix(`Smat')
	matname `Smat' naturalfreq Spectrum, c(.) e
	ret mat S = `Smat'
	mkmat `omega' `wave' `P', matrix(`Pmat')
	matname `Pmat' FourierFreq Period Periodogram, c(.) e
	ret mat P = `Pmat'
	if "`cinterval'" != ""{
		mkmat `freq' `SL' `SU', matrix(`CI')
		matname `CI' naturalfreq Lower Upper, c(.) e
		ret mat CI = `CI'
	}
}
else if "`log'" != ""{
	mkmat `freq' `S' `logS', matrix(`Smat')
	matname `Smat' naturalfreq Spectrum LogSpectrum, c(.) e
	ret mat S = `Smat'
	mkmat `omega' `wave' `P' `logP', matrix(`Pmat')
	matname `Pmat' FourierFreq Period Periodogram LogPergram, c(.) e
	ret mat P = `Pmat'
	if "`cinterval'" != ""{
		mkmat `freq' `SL' `SU', matrix(`CI')
		matname `CI' naturalfreq Lower Upper, c(.) e
		ret mat CI = `CI'
	}
}
matname `wvec' weights, c(.) e
ret mat W = `wvec'
ret scalar N = `n'
ret scalar dof = `dof'


if "`out'" != ""{
	keep `k' `freq' `omega' `wave' `Preal' `Pimaginary' `P' `S'
	rename `k' _k
	rename `freq' naturalfreq
	rename `omega' FourierFreq
	rename `wave' Period
	rename `Preal' Preal
	rename `Pimaginary' Pimaginary
	rename `P' Periodogram
	rename `S' Spectrum

	capture confirm file "`out'.dta"
		if _rc==0 {
			if "`replace'" != "" {
				qui save `"`out'"', replace
			}
			else {
				di as error "File `out'.dta already exists"
				exit
			}
		}
		else {
			qui save `"`out'"'
		}
}

restore

end

program define _convdaniell
	args c1 c2
	local m = 2*`c1'+1
	local n = 2*`c2'+1
	tempname a b c
	mat `a' = J(`m',1,1/`m')
	mat `b' = J(`n',1,1/`n')
	local h = max(`m'+`n'-1,`m',`n')
	mat c = J(`h',1,0)
	forv k=1(1)`h' {
		local sum = 0
		forv j=1(1)`m' {
			if `k'>=`j' & (`k'-`j'+1)<=`n' {
			local sum = `sum' + `a'[`j',1]*`b'[`k'-`j'+1,1]
			}
		}
		mat c[`k',1] = `sum'
	}

end

program define _convvector
	args u v
	local m = rowsof(`u')
	local n = rowsof(`v')
	local h = max(`m'+`n'-1,`m',`n')
	mat c = J(`h',1,0)
	forv k=1(1)`h' {
		local sum = 0
		forv j=1(1)`m' {
			if `k'>=`j' & (`k'-`j'+1)<=`n' {
			local sum = `sum' + `u'[`j',1]*`v'[`k'-`j'+1,1]
			}
		}
		mat c[`k',1] = `sum'
	}
end

version 12
mata:
void smoothdens(string scalar varname1, string scalar varname2)
	{
	 w = st_matrix(st_local("wvec"))
	 w = w/colsum(w)
	 q = rows(w)
	 b = (q-1)/2
	 Pmat = st_data(., varname1)
	 real matrix S
	 st_view(S=.,., varname2)
	 m = rows(Pmat)
	 P = Pmat[(1..m),.]
	 P = (P[(b+1..2),.]\P\P[(m-1..m-b),.])

	 for(i=1+b; i<=m+b; i++) {
		S[i-b,.] = w'*P[(i-b..i+b),.] :/(4*pi())
	 }
	 st_numscalar("r(wssq)",w'*w)
	 st_replacematrix(st_local("wvec"),w)
	 }
end