//
//  ProtoCD.ox
//  Sep 99 (PJW)
//
//  Ox version of the Cobb-Douglas PROTO model.  Solves for the base case
//  and one or more policy solutions (see the main routine).  Requires the
//  Newton's Method solution algorithm in newton.ox.
//

#include <oxstd.h>
#include "newton.ox"


//  Here are the parameters:

decl a     = 0.285;
decl alpha = 2    ;
decl g     = 0.5  ;


//  Here are the base case variables:

decl h  = 100.0;
decl p  = 1.000;
decl tc = 0.200;
decl tx = 0    ;
decl tw = 0    ;


//  These are the endogenous variables:

decl y, w, s, c, j, l, x, q, pc, px, wh, u, wl;
decl xlbl={ "Y", "W", "S", "C", "J", "L", "X", "Q", "Pc", "Px", "Wh", "U", "WL" };


//  These are used for computing EVs after the base run


decl base_pc;   // base case after-tax price faced by households
decl base_wh;   // base case after-tax wage received by households
decl base_y;    // base case full income


//  The following function evaluates the model for a given guess of
//  the endogenous variables in vector "guess".

f(guess)
{
   decl miss;

   //  Initialize the result vector so that it will have the
   //  correct number of elements.  Set each element to 1 to
   //  help detect mistakes in equation numbering below.

   miss = ones(13,1);

   //  Map the guess vector into the model's variables.

   y  = guess[ 0];
   w  = guess[ 1];
   s  = guess[ 2];
   c  = guess[ 3];
   j  = guess[ 4];
   l  = guess[ 5];
   x  = guess[ 6];
   q  = guess[ 7];
   pc = guess[ 8];
   px = guess[ 9];
   wh = guess[10];
   u  = guess[11];
   wl = guess[12];

   //  Evaluate the equations as LHS - RHS.  If we had accidentally
   //  left one out we'd be able to detect that because it would be
   //  set to 1 and the model would fail to solve.

   miss[ 0] = y    - wh*h - s                               ; 
   miss[ 1] = pc*c - a*y                                    ;
   miss[ 2] = wh*j - (1-a)*y                                ;
   miss[ 3] = wh*l - a*y + s                                ;
   miss[ 4] = l    - (q/alpha)*( (1-g)*px/(g*w) )^g         ; // modified
   miss[ 5] = x    - (q/alpha)*( g*w/((1-g)*px) )^(1-g)     ; // modified
   miss[ 6] = p*q  - w*l - px*x                             ; // modified
// miss[ 6] = p    - (1/alpha)*( px/g )^g*( w/(1-g) )^(1-g) ; // alternative
   miss[ 7] = s    - tc*p*c - tx*p*x - tw*w*l               ; // modified
   miss[ 8] = pc   - p*(1+tc)                               ; 
   miss[ 9] = px   - p*(1+tx)                               ; 
   miss[10] = wh   - w*(1-tw)                               ; 
   miss[11] = u    - (c^a)*(j^(1-a))                        ; // utility
   miss[12] = wl   - q + x + c                              ;

   return(miss);
}


//  A function to print out results nicely, including labels and a column of
//  percentage changes from the base case.

printResults(x,basex)
{
   decl column_formats = {"%9.5f","%10.4f %%"};
   decl ev, pct;
   decl labels;

   labels = xlbl;
   pct = 100*(x-basex)./basex;

   if( !(x == basex) )
      {
      ev     = u*(base_pc/a)^a*(base_wh/(1-a))^(1-a) - base_y ;
      x      = x | ev;
      pct    = pct | 100*ev/base_y;
      labels = labels | "EV";
      }

   print("%r", labels, "%cf", column_formats, x ~  pct );
}


main()
{
   decl x, basex;

   //  Initial guess of the endogenous variables.  This is a
   //  bad guess but not nearly as bad as zeros.

   x = ones(13,1);

   //  Experiment 1: Solve for the base case.

   println( "\nPROTO CD, Exp 1: Base Case\n" );
   x = solveNewton(f,x);
   printResults(x,x);

   //  Save the base case solution for use later.

   basex = x;

   base_pc = pc;
   base_wh = wh;
   base_y  = y;

   //  Experiment 2: Check that the model is homogeneous in wages and prices.

   p=1.5;

   println( "\nPROTO CD, Exp 2: Homogeneity Test, Expect 50s and 0s\n" );
   x = solveNewton(f,x);
   printResults(x,basex);

   p=1;

   //  Experiment 3: increase TC

   x=basex;
   tc=0.3;

   println( "\nPROTO CD, Exp 3: Policy Case, TC=0.3\n" );
   x = solveNewton(f,x);
   printResults(x,basex);

   tc=0.2;

   //  Experiment 4: increase TX

   x=basex;
   tx=0.1;

   println( "\nPROTO CD, Exp 4: Policy Case, TX=0.1\n" );
   x = solveNewton(f,x);
   printResults(x,basex);

   tx=0.0;

   //  Experiment 5: increase TW

   x=basex;
   tw=0.1;

   println( "\nPROTO CD, Exp 5: Policy Case, TW=0.1\n" );
   x = solveNewton(f,x);
   printResults(x,basex);

   tw=0.0;

}