optimize.c

//
// 
// 
// C functions for heavy duty
// 
// gcc -shared -fPIC -O3 optimize.c -lgsl -lm -lgslcblas -o findP.so
// icc -fast  -I/home/gcam/china/gsl -L/home/gcam/china/gsl/.libs -shared -fPIC optimize.c -lgsl -lgslcblas  -o findP.so 
// gcc -I/home/gcam/china/gsl -L/home/gcam/china/gsl/.libs -shared -fPIC optimize.c -lgsl  -o findP.so
// 

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_randist.h>

int linearIndex2d(int x, int nX, int y){ return x + nX*y; }
int linearIndex3d(int x, int nX, int y, int nY, int z){ return x + nX*(y + nY*z);}

int min(int x, int y){
	if(x>= y)
		return y;
	else
		return x;
}

double adjustmentCost(int kIn, int kPrimeIn, double delta, double parameterConvex, double parameterFixed,double* vK){
	// Perhaps non-convex costs

	// Fixed cost of adjustment
	double k = vK[kIn];
	double kPrime = vK[kPrimeIn];

	double changeInK = kIn == kPrimeIn ? 0 : 1; // to account for numerical imprecission
	return  changeInK*k*parameterFixed + ((k*parameterConvex)/2)*pow((kPrime-(1-delta)*k)/k,2); 

	// Capital Irreversibility
	// double disinvestment = kPrime > (1-delta)*k ? 0 : 1;
	// return disinvestment*(1-parameterFixed)*(kPrime - (1-delta)*k) + ((k*parameterConvex)/2)*pow((kPrime-(1-delta)*k)/k,2); 
}

double equityFloatationCosts(double dividends,  double parameter){

	return dividends > 0.0 ? 0.0 :  parameter*dividends;
}

double computeExpectationTruncatedLognormal(double cost, double mu, double sigma2){
	//http://home.datacomm.ch/paulsoderlind/Courses/OldCourses/EcmXSta.pdf
	double b0 = (log(cost) - mu)/sqrt(sigma2);
	return exp(mu + 0.5*sigma2)*gsl_cdf_ugaussian_P(-sqrt(sigma2)+b0)/gsl_cdf_ugaussian_P(b0);
}

struct costFinderParams {double depreciatedCapital; double netCashFlows; double* kGrid; int nK;  int nQ;
						 double continuation;  double floatationCost; double futureK; double expectedDebt;
						 double adjustmentCost;};

double difference(double cf, void* p){

	struct costFinderParams *params = (struct costFinderParams *)p;
	double dividends,costsOfFloatation,valueOfStaying,answer = 0.0;

	// Make sure the floatation cost takes into consideration the fixed cost of operation
	dividends = params->netCashFlows - cf;
	costsOfFloatation = equityFloatationCosts(dividends, params->floatationCost);
	// However, the value of staying is net of net production (as both exit and stay eat production net of debt repayment and wages)
	valueOfStaying = params->depreciatedCapital - params->futureK  + params->expectedDebt + params->continuation + costsOfFloatation - cf - params->adjustmentCost;
	answer = params->depreciatedCapital  - valueOfStaying;

	// printf("Guess cf %1.2f k(1-d) %1.2f divs %1.2f valOfStaying %1.2f  answer %1.2f flcost %1.5f \n",cf,params->depreciatedCapital,dividends,valueOfStaying,answer,equityFloatationCosts);

	return answer;
}

double findCutoffCost(struct costFinderParams *params ){
	
	double x_lo=0.001, x_hi = 5.0, answer = 0;
	int status;
	int iter = 0, max_iter = 1000;
	double root =  0;
	double bisectionTolerance = 1e-3;
	double minCost = 0.0; // smallest cf to consider for w - cf
	double maxCost = 1e5*params->kGrid[params->nK-1]; // largest cf to consider. Top of grid because if larger will for all entries violate 
										     // constraint that it has to be within the grid
			
	gsl_function F;
	F.function = &difference;
	F.params = params;

	const gsl_root_fsolver_type * T  = gsl_root_fsolver_bisection;
	gsl_root_fsolver * s   = gsl_root_fsolver_alloc (T);


	gsl_root_fsolver_set (s,&F,minCost,maxCost);					
	do{
		iter++;
		status = gsl_root_fsolver_iterate (s);
		root = gsl_root_fsolver_root (s);
		x_lo = gsl_root_fsolver_x_lower (s);
      	x_hi = gsl_root_fsolver_x_upper (s);
      	status = gsl_root_test_interval (x_lo, x_hi,
                                        bisectionTolerance, 0.00000);
	}
	while (status == GSL_CONTINUE && iter < max_iter);

	answer = root;
	gsl_root_fsolver_free (s);

	return answer;
}


void freePointer(double* p){
	free(p);
	
}




double* findPolicyInCPricesPreviousK(int kInitial,int qInitial, int z, double theta, double delta, double beta,double* vK,
					 int nK, double* vQ, int nQ,double* vZ, int nZ, double* currentValueZ, double* transition, double rBorrow,double fixedCostMean, double fixedCostVar,double capitalShare,
					  double floatationCost, double elasticityOfSubstitution, double aggOutput,double adjustmentCostParameter,int previousK,double adjustmentCostParameterFixed){

	int policyK = -1;
	int policyQ = -1;
	double probOfExit = 0; // probability of exit ( P(c_f >= c^*))


	int k,q,zIndex;
	int feasibleSolution = 0;
	int allConstrainedBefore = 0; // if all q choices are hitting the constraint at this value of k'
								  //  advance to the next k'
	double highestVal = -1000000000009;
	double newVal, continuation;
	double inDiffCost,expectedCost,costsOfFloatation; // cost that leaves firm indifferent between leaving or staying tomorrow if z=1
	double bestDividend; // dividend at optimum
	double kValue;
	
	double qValue; // value of state contingent debt
	double dividends;  // dividends
	double distributions; // takes into account floatation cost
	double bestCost; // for storing expected fixed cost
	double bestProbability; // for storing prob of exit
	double adjCost,expectedDebt;

	double kCameInToPeriod = vK[kInitial];
	double debtCameIntoPeriod = vQ[qInitial];

	double realOutput  = vZ[z]*pow(kCameInToPeriod,capitalShare);					
	double price       = pow(realOutput/aggOutput,-1/elasticityOfSubstitution);
	
	for(k = previousK; k < nK; k++){

		allConstrainedBefore = 0;
		kValue  = vK[k];
		adjCost = adjustmentCost(kInitial,k,delta,adjustmentCostParameter,adjustmentCostParameterFixed,vK);

		for(q = 0; q < nQ; q++){

			// these make sure at this guess of k', the collateral constraint  is not violated
			qValue = theta*(1-delta)*kValue < (1+rBorrow)*vQ[q] ? theta*(1-delta)*kValue : vQ[q]; 
			expectedDebt = qValue;
			// are we constrained for this k at all q?
			allConstrainedBefore = theta*(1-delta)*kValue < (1+rBorrow)*vQ[q] ? 1 : 0; 

			
			dividends = price*realOutput - (1+rBorrow)*debtCameIntoPeriod  + (1-delta)*kCameInToPeriod - 
						kValue + expectedDebt - adjCost;



			// First check if exits even at zero fixed cost
			costsOfFloatation =  equityFloatationCosts(dividends,floatationCost);
			if( (1-delta)*kCameInToPeriod > (1-delta)*kCameInToPeriod - kValue + expectedDebt + continuation + costsOfFloatation - adjCost){
				inDiffCost = -1; // if exit even at zero cost, make probability of exit equal to 1.0
			} else { // else solve the nonlinear equation
				struct costFinderParams params = {(1-delta)*kCameInToPeriod,dividends,vK,nK,nQ,continuation,
				floatationCost,kValue,expectedDebt,adjCost};
				// printf("Finding cost for k=%d k'=%d q = %d q' = %d dividends at zero %f\n",kInitial,k,qInitial,q,dividends);
				inDiffCost = findCutoffCost(&params);
			}

			// If stay at positive fixed cost, exit probability is the probability that cost is >= to indifference cost. 
			// Otherwise exit with certainty

			probOfExit = inDiffCost > 0 ? 1-gsl_cdf_lognormal_P(inDiffCost,fixedCostMean,sqrt(fixedCostVar)) : 1.0;
			expectedCost = inDiffCost > 0 ? computeExpectationTruncatedLognormal(inDiffCost,fixedCostMean,fixedCostVar) : -5;
			dividends =  dividends - expectedCost;
			distributions = dividends + equityFloatationCosts(dividends,floatationCost);

			continuation = 0;
			for(zIndex = 0; zIndex < nZ; zIndex++){
				continuation += transition[linearIndex2d(z,nZ,zIndex)]*currentValueZ[linearIndex3d(k,nK,q,nQ,zIndex)];
			}

			continuation = beta*continuation;

			newVal = probOfExit*(price*realOutput - (1+rBorrow)*debtCameIntoPeriod + (1-delta)*kCameInToPeriod)+(1-probOfExit)*(distributions + continuation);
			// printf("k %d q %d z%d, q' %f has value %f allC %d  LHS %f RHS %f \n",kInitial,qInitial,z,qValue,newVal,allConstrainedBefore,theta*(1-delta)*kValue , vQ[q]);
			if(newVal > highestVal){

				highestVal = newVal;
				policyK         = k + 1;         //Julia 1 based indexing
				policyQ         = q + 1;
				bestProbability = probOfExit;
				bestDividend    = distributions;
			}
			if(allConstrainedBefore > 0)
				break;
		}
	}

	// printf("WOOHOO OUT k %d q %d z %d  vf %1.2f pk %d pq %d dist %1.2f \n",kInitial,qInitial,z,highestVal,policyK,policyQ,bestDividend);

    double *results = (double *)malloc(sizeof(double)*6);
    results[0] = policyK;
    results[1] = policyQ;
    results[2] = highestVal;
    results[3] = bestProbability;
    results[4] = bestDividend;
    results[5] = bestCost;
    return results;
}



double* findPolicyInCPrices(int kInitial,int qInitial, int z, double theta, double delta, double beta,double* vK,
					 int nK, double* vQ, int nQ,double* vZ, int nZ, double* currentValueZ, double* transition, double rBorrow,double fixedCostMean, double fixedCostVar,double capitalShare,
					  double floatationCost, double elasticityOfSubstitution, double aggOutput,double adjustmentCostParameter,double adjustmentCostParameterFixed){

	int policyK = -1;
	int policyQ = -1;
	double probOfExit = 0; // probability of exit ( P(c_f >= c^*))


	int k,q,zIndex;
	int feasibleSolution = 0;
	int allConstrainedBefore = 0; // if all q choices are hitting the constraint at this value of k'
								  //  advance to the next k'
	double highestVal = -1000000000009;
	double newVal, continuation;
	double inDiffCost,expectedCost,costsOfFloatation; // cost that leaves firm indifferent between leaving or staying tomorrow if z=1
	double bestDividend; // dividend at optimum
	double kValue;
	
	double qValue; // value of state contingent debt
	double dividends;  // dividends
	double distributions; // takes into account floatation cost
	double bestCost; // for storing expected fixed cost
	double bestProbability; // for storing prob of exit
	double adjCost,expectedDebt;

	double kCameInToPeriod = vK[kInitial];
	double debtCameIntoPeriod = vQ[qInitial];

	double realOutput  = vZ[z]*pow(kCameInToPeriod,capitalShare);					
	double price       = pow(realOutput/aggOutput,-1/elasticityOfSubstitution);
	
	for(k = 0; k < nK; k++){

		allConstrainedBefore = 0;
		kValue  = vK[k];
		adjCost = adjustmentCost(kInitial,k,delta,adjustmentCostParameter,adjustmentCostParameterFixed,vK);

		for(q = 0; q < nQ; q++){

			// these make sure at this guess of k', the collateral constraint  is not violated
			qValue = theta*(1-delta)*kValue < (1+rBorrow)*vQ[q] ? theta*(1-delta)*kValue : vQ[q]; 
			expectedDebt = qValue;
			// are we constrained for this k at all q?
			allConstrainedBefore = theta*(1-delta)*kValue < (1+rBorrow)*vQ[q] ? 1 : 0; 
			dividends = price*realOutput - (1+rBorrow)*debtCameIntoPeriod  + (1-delta)*kCameInToPeriod - 
						kValue + expectedDebt - adjCost;




			// First check if exits even at zero fixed cost
			costsOfFloatation =  equityFloatationCosts(dividends,floatationCost);
			if( (1-delta)*kCameInToPeriod > (1-delta)*kCameInToPeriod - kValue + expectedDebt + continuation + costsOfFloatation - adjCost){
				inDiffCost = -1;
			} else { // else solve the nonlinear equation
				struct costFinderParams params = {(1-delta)*kCameInToPeriod,dividends,vK,nK,nQ,continuation,
				floatationCost,kValue,expectedDebt,adjCost};
				inDiffCost = findCutoffCost(&params);
			}

			// If stay at positive fixed cost, exit probability is the probability that cost is >= to indifference cost. 
			// Otherwise exit with certainty

			probOfExit = inDiffCost > 0 ? 1-gsl_cdf_lognormal_P(inDiffCost,fixedCostMean,sqrt(fixedCostVar)) : 1.0;
			expectedCost = inDiffCost > 0 ? computeExpectationTruncatedLognormal(inDiffCost,fixedCostMean,fixedCostVar) : -5;
			dividends =  dividends - expectedCost;
			distributions = dividends + equityFloatationCosts(dividends,floatationCost);


			continuation = 0;
			for(zIndex = 0; zIndex < nZ; zIndex++){
				continuation += transition[linearIndex2d(z,nZ,zIndex)]*currentValueZ[linearIndex3d(k,nK,q,nQ,zIndex)];
			}

			continuation = beta*continuation;

			newVal = probOfExit*(price*realOutput - (1+rBorrow)*debtCameIntoPeriod + (1-delta)*kCameInToPeriod)+(1-probOfExit)*(distributions + continuation);
	
			if(newVal > highestVal){

				highestVal = newVal;
				policyK   = k  + 1;         //Julia 1 based indexing
				policyQ   = q + 1;
				bestProbability = probOfExit;
				bestDividend = distributions;
				bestCost = expectedCost;
			}
			if(allConstrainedBefore > 0)
				break;
		}
	}

	 // printf("WOOHOO OUT k %d q %d z %d  vf %1.2f pk %d pq %d dist %1.2f \n",kInitial,qInitial,z,highestVal,policyK,policyQ,bestDividend);

    double *results = (double *)malloc(sizeof(double)*6);
    results[0] = policyK;
    results[1] = policyQ;
    results[2] = highestVal;
    results[3] = bestProbability;
    results[4] = bestDividend;
    results[5] = bestCost;
    return results;
}