*! glmcstar Esarey c* loss aversion JEsarey NDanneman 10sep2009
program glmcstar
	version 10.1
	args Gamma

if `: word count `0'' != 1 {					//displays an error if user inputs any number of arguments other than one
mata _error("cstar takes exactly one integer as an argument")
}			

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

quietly glmcstargears `Gamma'

local NumBetas = colsof(CSTAR)
local regs : colnames e(b)


display as text "c* values for the previous regression"	//displays C-star generic info			
display as text "gamma = `Gamma'"				//reports the gamma parameter input by the user

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

forvalues i = 1/`NumBetas' {
display as text word("`regs'", `i') 				//displays the name of each regressor in green
display as result CSTAR[1,`i']					//displays each c-star value in yellow
di ""								//adds a blank line between a c-star value and the following regressor's name
}

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

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

end

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

program glmcstargears
	version 10.1
	args Gamma

matrix b = e(b)							//creates matrix b, a rowvector of betas
matrix V = e(V)							//V, a rowvector of variances


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

end

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

mata:

void function CstarExecute()
{
Betas = st_matrix("b")					//pulls the beta matrix in from stata
AbsBetas = abs(Betas)					//absolute values of beta matrix
Sign = sign(Betas)					//sign of the beta matrix
Variances = diagonal(st_matrix("V"))			//same for variances, but it is a column
Variances = Variances'					//makes Variances a rowvector
Sigmas = sqrt(Variances)	                      	//turns Variances into Sigmas via square root
Pieces = (AbsBetas \ Sigmas \ Sign)			

NumberOfBetas = cols(Pieces)				//to get the number of times to loop
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 rowvector Output					//declare Output
Output = J(1,NumberOfBetas,.)				//set up empty Output rowvector
for (i=1; i<=NumberOfBetas; i++) {
	m = Pieces[1,i]
	sigma = Pieces[2,i]
	sign = Pieces[3,i]
	lower = m - (8*sigma)				//extended the limits to 8*sigma from 6
	higher = m + (8*sigma)
	Output[1,i] = sign * rootfinder(lower, higher, numslic, a, m, sigma)			//fill empty Output vector with results of rootfinder
	}

cstarsign = sign(Output)				//creats a vector of the signs of the output

for (i=1; i<=NumberOfBetas; i++) {				//this loop assigns c* values of zero to those that would otherwise have the opposite sign of their betas
	if (cstarsign[1,i] != Pieces[3,i]) Output[1,i] = 0
	}

st_matrix("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