//
//  Proto.ox
//  Sep 99 (PJW)
//
//  Ox version of the 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 b = 2.040;


//  Here are the base case variables:

decl h = 100.0;
decl t = 0.200;
decl p = 1.000;


//  These are the endogenous variables:

decl y, w, s, c, j, l, q, wl;
decl xlbl={ "Y", "W", "S", "C", "J", "L", "Q", "WL" };

//  The following function evaluates the model:

f(x)
{
   decl val;

   //  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.

   val = ones(8,1);

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

   y  = x[0];
   w  = x[1];
   s  = x[2];
   c  = x[3];
   j  = x[4];
   l  = x[5];
   q  = x[6];
   wl = x[7];

   //  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.

   val[0] = y - w*h - s      ;
   val[1] = p*(1+t)*c - a*y  ;
   val[2] = w*j - (1-a)*y    ;
   val[3] = w*l - a*y + s    ;
   val[4] = l - q/b          ;
   val[5] = p - w/(b-1)      ;
   val[6] = t*p*c - s        ;
   val[7] = wl - q + q/b + c ;

   return(val);
}


//  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 = {"%7.3f","%8.2f %%"};

   print("%r", xlbl, "%cf", column_formats, x ~  100*(x-basex)./basex );
}


//  The main routine

main()
{
   decl x, basex;

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

   x = ones(8,1);

   //  Solve for the base case.

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

   //  Save the base case solution for use later.

   basex = x;

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

   p=2;

   println( "\nPROTO, Homogeneity Test, Expect 100s and 0s\n" );
   x = solveNewton(f,x);
   printResults(x,basex);

   p=1;

   //  Check that the model is homogeneous in quantities

   h=200;

   println( "\nPROTO, Quantity Homogeneity Test, Expect 100s and 0s\n" );
   x = solveNewton(f,x);
   printResults(x,basex);

   h=100;

   //  Now change the tax rate and solve for the new equilibrium.
   //  Since we didn't change X we will begin from the old solution,
   //  which is generally a good idea.

   x=basex;
   t=0.3;

   println( "\nPROTO, Policy Case, T=0.3\n" );
   x = solveNewton(f,x);
   printResults(x,basex);

   t=0.2;
}