*! glmcstari cstar c* loss aversion Esarey JEsarey NDanneman 10sep2009
program glmcstari
	version 10.1
	args Beta SE Gamma

if `: word count `0'' != 3 {					//displays an error if user inputs any number of arguments other than three
mata _error("glmcstari takes exactly 3 integers as arguments")
}			

if `Gamma' ==0 {						//displays an error if user inputs 0 as the gamma parameter
mata _error("the loss aversion parameter cannot equal zero")
}			

quietly glmcstarigears `Beta' `SE' `Gamma'				//run the guts of the program

display ""							//adds a blank line for clarity
display as text "c* value for the conditions:"			//display this as a header						//	***
display ""		
display as text "gamma = `Gamma'"				//reports the gamma parameter input by the user
display as text "Beta = `Beta'"
display as text "SE = `SE'"

display ""							//adds a blank line to reduce visual clutter

display as result CSTAR						//displays the c* value in yellow
di ""								//adds a blank line to reduce visual clutter


if `Gamma' == 1 {
display "Note: you have entered a loss aversion factor of 1, so cstar has returned the original beta."		
}

if `Gamma' < 1 {
display "Note: you have entered a loss aversion factor less than one, signifying risk-seeking behavior."	
}

end

**********************************************

program glmcstarigears
	version 10.1
	args Beta SE Gamma

mata: CstarExecute()						//executes the mata function CstarExecute

end

*********************************************

mata:

void function CstarExecute()
{
BetaStr = st_local("Beta")				//pulls the beta in from stata
BetaM = strtoreal(BetaStr)				//makes the beta a real instead of a string
SEStr = st_local("SE")					//pulls the SE in from stata
SEM = strtoreal(SEStr)					//makes the SE into a real (from a string)
AbsBeta = abs(BetaM)					//absolute value of beta
Sign = sign(BetaM)					//sign of beta


numslic = 5000						//extended this for higher accuracy
Gammaparm = st_local("Gamma")				//this is treated as a string scalar, but the program needs a real scalar
a = strtoreal(Gammaparm)				//convert string to real, and give it correct letter

real scalar Output					//declare Output

m = AbsBeta
sigma = SEM
sign = Sign
lower = m - (8*sigma)					//extended the limits to 8*sigma from 6
higher = m + (8*sigma)

Output = sign * rootfinder(lower, higher, numslic, a, m, sigma)				//fill empty Output scalar with results of rootfinder
	
cstarsign = sign(Output)
if (Sign != cstarsign) Output = 0			//this command sets c* equal to zero if it has the opposite sign of its beta


st_numscalar("CSTAR", Output)

}

//*********************************
//THE UTILITY FUNCTION

function utility(
real scalar x,								//variable of integration
real scalar a,									//gamma - the loss parameter
real scalar m,								//beta
real scalar sigma,							//sd of beta
real scalar t)								//variable, roots taken over t
{

k = ( (x>t) ? 1 : -1 )							//a logical expression

distr = ( normalden(x, m, sigma) )

return( (a^(-k)) * (distr) * (x-t) )

}

//*********************************
//INTEGRATION
//An n-slice integration routine using simpson's approximation, use 120 slices for numslic

scalar NDint(	                                                    	     	 	// to perform integration via simpson'sapproximation
real scalar t,
pointer scalar f,
real scalar lower,                                                                      //lower bound of integration
real scalar higher,                                                                     //upper bound
real scalar numslic,									// where numslic = number of slices, must be an EVEN number
real scalar a,
real scalar m,
real scalar sigma)
{

 if (lower > higher) _error("limits of integration reversed")        			//make sure upper > lower NB: better to use _error() than return()
 if ( (numslic/2) != (floor(numslic/2)) ) _error("numslic must be an EVEN number")	//ensures that the number of slices used is even
 
 real scalar width									//declarations, for assurance's sake
 real colvector Y
 real rowvector X
 real scalar numslicplsone								//slices plus one
 numslicplsone = numslic +1

 width = (higher-lower) / numslic                                                       //calculates width of each "slice"

 Y = J(numslicplsone,1,.)                                                               //creates an empty column vector
 for (i=1; i<=numslicplsone; i++) {
    x = lower + ((i-1)*width)
    Y[i,1]  = (*f)(x,a, m, sigma, t)                                               	//places function evaluations directly into Y
 }


 X = J(1,numslicplsone,.)						//this loop creates a vector of the fudge factors for integration via simpson's rule
 for (i=1;i<=numslicplsone; i++) {					// ie. (1, 4, 2...2, 4, 1)
  if (i==1) X[1,i] = 1
  if (i==numslicplsone) X[1,i] = 1  
  if ((i/2) == (floor(i/2))) X[1,i] = 4
  if ( ((i+1)/2) == (floor( ((i+1)/2)) ) & (i!=1) & (i!=numslicplsone) ) X[1,i] = 2
}

return(X*Y*width/3)							  //matrix multiplication makes this quicker and more readable
}

//*********************************
//THE ROOT FINDER - uses mm_root from the moremata package

function rootfinder(
real scalar lower,
real scalar higher,
real scalar numslic,
real scalar a,
real scalar m,
real scalar sigma)
{

mm_root(answer=., &NDint(), -m, m, 0, 1000, &utility(), lower, higher, numslic, a, m, sigma)

return(answer)

}

end