/*This ado file generates random data according to the sgt distribution.

Author --Jacob Orchard

v1.0
Update 7/29/2016*/




program rsgt, eclass
version 13.0
	if replay() {
		display "Replay not implemented"
	}
	else {
	args varname mu lambda sigma p q 
		confirm name `varname'
		confirm number `mu' 
		confirm number `lambda'
		confirm number `sigma'
		confirm number `p'
		confirm number `q'
		
		if `p' <=0{
			di as error "Parameter p must be positive"
		}
		if `q' <=0{
			di as error "Parameter q must be positive"
		}
		if `sigma' <= 0{
			di as error "Parameter sigma must be positive"
		}
		if `lambda' <= -1 | `lambda' >= 1{
			di as error "Parameter lambda must be between -1 and 1."
		}
		
		quietly mata: rsgt("`varname'",`mu',`lambda',`sigma',`p',`q')
	}

end



//Mata functions

version 13
mata:
	function rsgt(string myvar, scalar mu, scalar lambda, scalar sigma, scalar p, scalar q)
	{
	nobs = st_nobs()
	base = runiform(nobs,1)
	newvar = J(nobs,1,0)
	paravec = mu, lambda, sigma,p,q
	
	if (nobs <1000){
	
		for (i=1; i <=nobs; i++){
			newvar[i] = estimatesgt(paravec,base[i])
			}
		}
	else{
		cdflist = interpolatesgt(paravec)
		
		for (i=1; i <=nobs; i++){
			newvar[i] = nearest_sgt(cdflist,base[i])
			}
	
	}
	
	(void) st_addvar("double", myvar)
	st_store(.,myvar,newvar[.,1])
	}
	end
	
mata:
	function estimatesgt(paravec,unifval)
	{
	
		mu      = paravec[1]
		lambda  = paravec[2]
		sigma   = paravec[3]
		p       = paravec[4]
		q       = paravec[5]
		
		S = optimize_init()
		optimize_init_which(S,"min")
		optimize_init_technique(S,"dfp nr")
		optimize_init_evaluator(S,&closestsgt())
		
		//Expected Value (starting value)
		start = mu + (sigma/2)*((q^(1/p)*exp(lngamma(1/p) + lngamma(q-(1/p)) - lngamma(q)))/exp(lngamma(1/p) + lngamma(q-(1/p)) - lngamma(q)))*((1+lambda)^2 -(1-lambda)^2)
		optimize_init_params(S,start)
		optimize_init_argument(S,1,paravec)
		optimize_init_argument(S,2,unifval)
		_optimize(S)
	    ans    = optimize_result_params(S)
		return(ans)
	
	}
	
end

mata:
	function interpolatesgt(paravec)
	{
		cdflist = J(1000,1,0)
		
		for (i=1; i <=1000; i++){
		cdflist[i] = estimatesgt(paravec,(i/1000))
		}
	return(cdflist)
    }
end

mata:
	function nearest_sgt(cdflist,unifval)
	{
		nx = 1000*unifval
		intx = floor(nx)
		remx = nx - intx
		next = intx+1
		
		if (intx ==0) {
			closesgt = cdflist[1]*remx
		}
		else{
			closesgt = cdflist[intx]*(1-remx) + cdflist[next]*remx
		}
		return(closesgt)
	}
end

mata: 
	function closestsgt(real scalar todo, real rowvector p, paravec,unifval, v,g,H)
	{
	
	v = abs(sgt_cdf(p,paravec)-unifval)
	
	}
	end
	
mata:
function sgt_cdf(matrix y, paravec)
	{
	ones = J(1,cols(y),1)
	mu     = paravec[1]
	lambda = paravec[2]
	sigma  = paravec[3]
	p      = paravec[4]
	q      = paravec[5]
	
	
	
	zu = abs(y-mu):^p/(abs(y - mu):^p + q*sigma:^p*(1+lambda*sign(y-mu))^p)
	F = .5*(1-lambda) + .5*(1+lambda*sign(y-mu))*sign(y-mu)*exp(lngamma(1/p) + lngamma(q) - lngamma(1/p+q))*ibeta(1/p,q,zu)
	return(F)
	}
end