/*
  ExpectedCV.ox
 Version history
 1. Anders Karlstrom November 2000
 2. Anders Karlstrom and Edward Morey, updated April 19. 2003. This version does not use differentiation of choice probs.


 The program calculates the expected cv using the "expected expenditure formula (Karlstrom and Morey (2003))using the estimated
 discrete-choice random utility model described in Karlstrom and Morey (2003), hereafter referred to as the "Maine model".

 More details on the application can be found by going to http://www.colorado.edu/Economics/morey/dataset.html. It lists the
 articles that use the "Maine" data set

 Elsewhere, the data was used to estimate the maximum likelihood estimates
 of the parameters in that model

 This program reads in those parameter estimates and the data that generated them.

 The output from this program consists of the vector of expected cv's for each individual in the sample and the sample average for
 the specified change in prices and quality
 

 The reader will need to modify this code for their specific application

 (users of gauss should find the code moderately familiar - Ox can easily read in gauss data sets, matrices and even Gauss code)

 Modifications required for other application include, but are not limited to,
 reading in the appropriate parameter estimates and data,
 correctly specifying the number of alternatives in the choice set,
 specifiying the conditional indirect utility functions,
 and correctly specifiying the change to be valued.

 This code works for both models with and without income effects (if no income effects are assumed, the results are the logsum
 results)


 This program creates a number of "functions". Function are commands that have inputs and output. For example, the built-in
 Ox command "log()" calculates the natural log.

 
 This program calls a number of functions that were defined in other files (Romberg_.ox and Solve_.ox)
 These include, for example, the created functions "zbrac" and "rtbis".

 In terms of the code below, most of it creates, section by section, a number of functions.

 The actual main program, "main()", to calcuale the expected cv's is quite short and is found at the bottom of this file.
 */

#include <oxstd.h>		  // this command appears in all Ox programs. It provides access to the Ox built-in functions


abs(const x)				// "abs" is a created function and returns the absolute value
{
return sqrt(sqr(x));
}

// The following functions are defined in the file Solve_.ox
zbrac(const tmp0,const tmp1, const retval,const zfunc);
rtbis(const mGuess, const t, const zfunc);
zfunc(const x);

/* for those unfamiliar with Ox, the term "const tmp0" simply says the first input into the function is a constant with the
temporary name tmp0. So, for example, zbrac has 4 elements (inputs plus outputs) */


/* The following variables and parameters are used by (inputs)into the functions defined in
Romberg_.ox and Solve_.ox  */
decl s_static=0;
decl ZFACTOR=1.6;
decl MAXNTRY=50;//50;
decl BISERR = 0.000001;
decl MAXBIS = 40;
decl rounderror = 1E-3;  // Sekant method is not exact, allow a little round off error, size dependent on size of income (<1)
decl ROMBEPS=0.1;

decl x,p0,q0,y0,p1,q1;  
decl b; 					   // b is the estimated parameter vector
decl mNOBS, mNrAlternatives;  // number of observations and number of alternatives
decl my, mV0;
decl mI ;					  // This is a global variable sent from Integrate() to integrand()
decl tid;					  // Used by the clock to measure elapsed time

#include "Romberg_.ox"	 // We need to include all the functions defined in Romberg_.ox



//-------------------------------------------------------------------------------------------------
// Read in your data here.
// Attributes are (p,q). Original attributes (prices and qualities) are read into p0 and q0 
// Original incomes should be read into y0
// Note that y0, p0 and q0 are used by the program elsewhere, so these need to be defined like this
ReadData()
{
  x=loadmat("x.mat");		// Data (not parameter estimates) is in a file named "x" with the .mat format - see comment below
  p0=x[][14:21];  p1=p0;	// Original prices p0. In our case columns 14-21 of the x matrix. Initially set p1=p0
  q0=x[][6:13];	  q1=q0;	// Original quality attributes q0. In our case, columns 6-13 of the x matrix
  y0 = x[][4] / 50;			// Original income y0 is in colum 4 (it is divided by 50 to express it per each of the assumed 50 choice occasions
  mNOBS=sizer(x);			// number of observations is the number of rows in x
  mNrAlternatives=9;		// Set the number of alternatives, We have 9: 8 fishing sites and nonparticipation

/*
explaining the above notation "x[][]" indicates which rows and columns of x to capture. For example,
x[][14:21] captures all of the rows but only columns 14 thru 21

Our application has 8 fishing alternative and nonparticipation, so 8 prices (columns 14:21)are read in.
Each fishing site has only one quality variable, catch. The 8 catch rates (columns 6:13)are read in for
each angler

Note that the type of file read in by the function "loadmat"
depends on the extension of the file name: 
.mat: ASCII matrix file, described below, 
.dat: ASCII data file with load information, 
.in7: PcGive 7 data file (with corresponding .bn7 file), 
.xls: Excel worksheet or workbook file, 
.wks and .wk1: Lotus spread sheet file, 
.dta: Stata data file (version 4--6), 
.dht: Gauss data file (v89, with corresponding .dat file), 
.fmt: Gauss matrix file (v89 and v96); 
A matrix file holds a matrix, preceded by two integers which specify the number
of rows and columns of the matrix.It will normally have the .mat extension.

The estimated parameter vector is read in by the function "main()", defined below

For those familiar with Gauss, note that Gauss starts indexing vectors by 1, whereas Ox (as C++) starts with 0.

*/
  
}

//-------------------------------------------------------------------------------------------------


/*
This section is where the user specifices the policy change(p1,q1) that is to be valued.

It defines a function called "PolicyChange"

For example, this program, as set up, halves the catch rate at the
first site. Other possible changes are coded below, but commented out.


PolicyChange() will have to be modified for other applications
*/

PolicyChange()
{							  

  q1[][0]= .5 * q0[][0];		// Ex: only one change. Halving catch rate of first alternative.
								// Again note that in Ox we index first column with 0.
/*
notationally, [0] refers to the first element. E.g. q0[][0] is the first column of the matrix q0

To value a different scenario, replace "q1[][0]= .5 * q0[][0];" with, for example, one of the following
scenarios

	q1[][0]= 1.1* q0[][0];   	  	 // Ex: increase catch rate of two first sites by 10%
    q1[][1]= 1.1* q0[][1];	    
								
    q1[][0]= 1.1* q0[][0];	  	  // Ex: increase catch rate of first site by 10%, and decreases it by 50% at second site
    q1[][1]= .5* q0[][1];	  
	
    p1[][0]= 1.1* p0[][0];	// Ex: Increase prices on three first alternatives by 10%  
    p1[][1]= 1.1* p0[][1];	  
    p1[][2]= 1.1* p0[][2];	  
*/

}

//-------------------------------------------------------------------------------------------------

/* The program needs to input the functional forms that specify the deterministic part of the
conditional indirect utility function for each alternative.

A function "Indirects()" is defined with three inputs: price, catch rate and income.

This section will need to be modified by the user

For ease of specification, the income argument is specified as a vector, one for each alternative
(hence, the income argument may be an NxJ matrix, where N is nobs and J is the number of alternatives.
So, each element in a row of the income matrix is identical.									  

The integration routine does not use this matrix feature of the income argument, but Root() does.

in: (p,q,y), NxM matrices with prices, qualities and incomes
out: NxM matrix with indirect utilities V(p,q,y)
			where N is nobs and J is the number of alternatives.

Indirects() will need to be modified by the user
*/


Indirects(const s_p, const s_q, const m_y)

{

// the dec1 lines name the local variables used by the function Indirects
  decl j, mV, b11;
  decl mV1, mV2, mV3, mV4, mV5,mV6,mV7,mV8,mV0;


  decl pp,pm,pd,pk,ps,pns,pnb,pq,cp,cm,cd,ck,cs,cns,cnb,cq, ppy;
  decl y, m_n, years, age, club;
  decl tp, tm, td,tk,ts,tns, tnb, tq;
  
  pp=x[][14];	pm=x[][15]; pd=x[][16];pk=x[][17]; ps=x[][18]; pns=x[][19]; pnb=x[][20]; pq=x[][21];
  // there are 8 fishing sites. For example pp is the price for the Penobscot river, and pq is for rivers in Quebec

  cp=x[][6];	cm=x[][7] ; cd=x[][8]; ck=x[][9] ; cs=x[][10]; cns=x[][11]; cnb=x[][12]; cq=x[][13];
  // cp is the catch rate at the Penobscot river
  tp=x[][22];	tm=x[][23]; td=x[][24];tk=x[][25]; ts=x[][26]; tns=x[][27]; tnb=x[][28]; tq=x[][29];
  // tp is, or example, the number of trips to the Penobscot

  years=x[][1]; club=x[][2]; age=x[][3]; y=x[][4]; m_n=x[][5];
  // years, for example, is the number of years of fishing experience


  b11=b[11];

  if (columns(m_y)>1)
  {
    ppy=m_y; 
  }
  else ppy=m_y .* ones(mNOBS,mNrAlternatives);

  // the following are the deterministic part of the conditional indirect utility functions for our 9 alteranatives

  mV0  = (b11 * (b[1][] * (ppy[][0] - s_p[][0]) + b[2][] * s_q[][0] + b[3][] * s_q[][0]  .^ .5 + b[10][] * ((1728.4720 + ppy[][0] - s_p[][0]) .^ .5)));    
  mV1  = (b11 * (b[1][] * (ppy[][1] - s_p[][1]) + b[2][] * s_q[][1] + b[3][] * s_q[][1] .^ .5 + b[10][] * ((1728.4720 + ppy[][1] -s_p[][1]) .^ .5)));    						;
  mV2  = (b11 * (b[1][] * (ppy[][2] - s_p[][2]) + b[2][] * s_q[][2] + b[3][] * s_q[][2].^ .5 + b[10][] * ((1728.4720 + ppy[][2] -s_p[][2]) .^ .5)));    
  mV3  = (b11 * (b[1][] * (ppy[][3] - s_p[][3]) + b[2][] * s_q[][3] + b[3][] * s_q[][3].^ .5 + b[10][] * ((1728.4720 + ppy[][3] -s_p[][3]) .^ .5)));    
  mV4  = (b11 * (b[1][] * (ppy[][4] - s_p[][4]) + b[2][] * s_q[][4] + b[3][] * s_q[][4].^ .5 + b[10][] * ((1728.4720 + ppy[][4] - s_p[][4]) .^ .5)));    
  mV5  = (b11 * (b[1][] * (ppy[][5] - s_p[][5]) + b[2][] * s_q[][5] + b[3][] * s_q[][5].^ .5 + b[10][] * ((1728.4720 + ppy[][5] - s_p[][5]) .^ .5)));    
  mV6  = (b11 * (b[1][] * (ppy[][6] - s_p[][6]) + b[2][] * s_q[][6] + b[3][] * s_q[][6].^ .5 + b[10][] * ((1728.4720 + ppy[][6] - s_p[][6]) .^ .5)));    
  mV7  = (b11 * (b[1][] * (ppy[][7] - s_p[][7]) + b[2][] * s_q[][7] + b[3][] * s_q[][7].^ .5 + b[10][] * ((1728.4720 + ppy[][7] - s_p[][7]) .^ .5)));    
  /* the next one is for nonparticipation */
  mV8 = (b[1][] * ppy[][8] + b[4][] + b[5][] * years + b[6][] * club + b[7][] * age + b[8][] * years .^ .5 + b[9][] * age .^ .5 + b[10][] * ((1728.4720 + ppy[][8]) .^ .5));
  
  mV= mV0 ~ mV1 ~ mV2 ~ mV3 ~ mV4 ~ mV5 ~ mV6 ~ mV7 ~ mV8;
  return mV;
}


//-------------------------------------------------------------------------------------------------

/* 
This section creates a function, "prob()" to calculate the predicted probabilities associated with each alternative

in: is the NxJ matrix, V, of the conditional indirects
out: is a NxJ matrix of choice probabilities

These choice probabilities correspond to our 2-level nested logit model.

prob() will need to be modified by the user to fit their discrete choice model (e.g nesting specification)
*/

prob(const V)

{
	decl mV, inclusm, inclusc, inclusp, inclus, linclus, lsump, lmsum, lcsum, ex, prob;

	mV=V;
	ex=exp(mV);			// input to this function: mV is indirect utilities
		
	inclusm = (sumr(ex[][0:4])) .^ (b[12][] / b[11][]);		// this is your code, basically
    inclusc = (sumr(ex[][5:7])) .^ (b[12][] / b[11][]);
	inclusp = (inclusm + inclusc) .^ (1 / b[12][]);
    /* inclus is the denomin in all the prob */
    inclus = ex[][8] + inclusp;		   
	linclus = log(inclus);
    lsump = ((1 / b[12][]) - 1) * log(inclusm + inclusc);
    lmsum = ((b[12][] / b[11][]) - 1) * log(sumr(ex[][0:4]));
    lcsum = ((b[12][] / b[11][]) - 1) * log(sumr(ex[][5:7]));

	prob= (mV[][0:4] + lsump +lmsum - linclus) ~ (mV[][5:7] + lsump +lcsum - linclus) ~ (mV[][8] - linclus);

	return exp(prob[][]);	// return choice probabilities
	
}


//-------------------------------------------------------------------------------------------------


/* This section creates the function, "Roots()"

It finds integral limits (solves for the roots).
That is, it calculates the needed income y_{ii} an individual needs when choosing alternative i both before and after the policy change.
It is a deterministic value, and is found by sloving V_i(p0,q0,y0) = V_i(p1,q1,y_{ii}) for each individual and each alternative


The output is returned by the global variable my, and is used in the function Integrate()

Roots() does not need to be modified for other application
*/

Roots()
{
  decl retval, mBrac,ytmp,tmp1,tmp0,j;
  mBrac = zbrac(0.2*y0,5*y0, &retval,zfunc);
  retval = rtbis(mBrac, &tmp1, zfunc);
  println("");
  if (retval==1) println("Found roots, integral limits ok");
  my=tmp1;

   if (isnan(my) || retval==0) {print("Generate: Error 546 Could not evaluate conditional compensations my_ii"); exit(9);}
}

//-------------------------------------------------------------------------------------------------


/* This section creates the function, "integrand()"

in: income at the compensated state
out: the survivial choice probability of alternative mI (a global variable)

integrand does not need to be modified for other applications
*/

integrand(const y)
{
  decl mV1, mVi, integrand;
  mV1=Indirects(p1,q1,y);
  mVi= mV0 .* (mV0 .>= mV1) + mV1 .* (mV1 .> mV0);
  mVi[][mI] = mV0[][mI];
  integrand =  prob(mVi)[][mI];

  return integrand;
}

//-------------------------------------------------------------------------------------------------

/* This section creates the function, "Integrate()"

It will do the integration to calculate the expected cv's

No inputs
out: the average of the expected compensating variation in the sample = 1/n sum_i E(cv_i)

Note that it calls the integration function "qromb" that is defined in Romberg_.ox 

Integrate() does not need to be modified for other applications
*/

Integrate()
{
  decl tmp,a,aa,prb,values, CVi, i;
  
  CVi=zeros(mNOBS,1);
  mV0=Indirects(p0,q0,y0);	// Original utilities

   for (i = 0; i < mNrAlternatives; ++i)      
  {
	values=zeros(mNOBS,1);
	mI=i;				   // mI is a global variable that is used in the integrand(y) function. 
	aa=limits(my');		   
	a=aa[0][]';             // (a,aa) is lower and upper integral limits.
							// a and aa are vectors with number of rows equal to number of observations (nobs)
	aa=my[][i];
	if (abs(a-aa) <= 0.00001) values=0;	 // if zero integral interval, integral evaluates to zero
	else						
	  qromb(integrand,a,aa,&values, ROMBEPS);	// qromb() integrates a vector of functions (integrand) from a to b, and returns the value in vector values
	tmp=prob(mV0)[][mI] .*  a + values;
	CVi += tmp;	  
  }



  print("\n CVi", CVi - y0 );
  return meanc(CVi-y0);
}

#include "solve_.ox"

//-------------------------------------------------------------------------------------------------

/*
This, the final code, creates the main program, "main()"

The output from this main program is the estimated expected cv's using the expected expenditure function
It will also print out the the time it took do the calculations.

Note that this main program, main(), is very short because all of the work is being done by the functions created above.

*/

main()
{
  decl  res,tid;
  tid=timer();					// Start the clock

  /* This program needs to read the parameter vector for the estimated model.
  The vector "b_IncE.mat" provided with this program is the parameter vector with income effects
  To calculate E(cv) for another model with another parameter vector, the user should read in another parameter vector here. */

  b=loadmat("b_IncE.mat"); 	    // Parameter vector with income effects	See Output_1.txt

  ReadData(); 
  PolicyChange();
  Roots();
  res=Integrate();
  println("Elapsed time:",timespan(tid), " seconds");
  println("E(cv)=",res);
}