//----------------------------------------------------------------------//
//  Market.ox
//  Sep 99 (PJW)
//
//  Simple partial equilibrium model of multiple agents interacting
//  in a single market.  Demonstrates the use of solveNewton.ox
//
//    Consumer type i demand:  qd(i) = a(i)*pb^(-e(i))
//    Firm type i supply:      qs(i) = b(i)*ps^(h(i))
//----------------------------------------------------------------------//

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

//
//  Here are the parameters for different consumers and firms
//

decl a1  = 1.0;  // consumer type 1's demand scale parameter
decl e1  = 1.1;  // consumer type 1's demand elasticity

decl a2  = 2.0;  // consumer type 2's demand scale parameter
decl e2  = 1.5;  // consumer type 2's demand elasticity

decl b1  = 1.0;  // firm type 1's supply scale parameter
decl h1  = 1.0;  // firm type 1's supply elasticity

decl b2  = 3.0;  // firm type 2's supply scale parameter
decl h2  = 2.0;  // firm type 2's supply elasticity

//
//  Populations of each type of consumer and firm
//

decl nd1 = 100;  // number of type 1 consumers
decl nd2 =  50;  // number of type 2 consumers
decl ns1 = 100;  // number of type 1 firms
decl ns2 = 100;  // number of type 2 firms

//
//  Here is the base case exogenous variable
//

decl tax = 0.200;


//  These are the endogenous variables:

decl qd1, qd2, qs1, qs2, qd, qs, pb, ps;
decl xlbl={ "Qd1", "Qd2", "Qs1", "Qs2", "Qd", "Qs", "Pb", "Ps" };


//----------------------------------------------------------------------//
//  f()
//
//  This function evaluates the "miss distances" for a given guess of
//  the model's endogenous variables.
//----------------------------------------------------------------------//

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.

   qd1 = x[0];
   qd2 = x[1];
   qs1 = x[2];
   qs2 = x[3];
   qd  = x[4];
   qs  = x[5];
   pb  = x[6];
   ps  = 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] = qd1 - a1*(pb)^(-e1) ;
   val[1] = qd2 - a2*(pb)^(-e2) ;
   val[2] = qs1 - b1*(ps)^(h1)  ;
   val[3] = qs2 - b2*(ps)^(h2)  ;
   val[4] = qd - ( nd1*qd1 + nd2*qd2 );
   val[5] = qs - ( ns1*qs1 + ns2*qs2 );
   val[6] = qd - qs ;
   val[7] = pb - ps*(1+tax);

   return(val);
}


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

printResults(x,basex)
{
   decl percent, column_formats, output;

   column_formats = {"%7.3f","%8.2f %%"};

   percent = 100*(x-basex)./basex ;
   output  = x ~ percent;

   print("%r", xlbl, "%cf", column_formats, output );
}


//----------------------------------------------------------------------//
//  main()
//
//  This is the program itself.  It sets up the experiments and calls
//  solveNewton to find the solutions.
//----------------------------------------------------------------------//

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( "\nMARKET, Base Case\n" );
   x = solveNewton(f,x);
   printResults(x,x);

   //  Save the base case solution for later comparisons.

   basex = x;

   //  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;
   tax=0.3;

   println( "\nMARKET, Tax Increase to 30%\n" );
   x = solveNewton(f,x);
   printResults(x,basex);

   tax=0.2;

   //  See what happens if more type 1 firms appear.  Change ns1 and
   //  solve for a new equilibrium.  Start from the original x.

   x=basex;
   ns1=200;

   println( "\nMARKET, Doubling Number of Type 1 Firms\n" );
   x = solveNewton(f,x);
   printResults(x,basex);

   ns1=100;


}