//----------------------------------------------------------------------//
//  harberger4.ox
//  Mar 04 (PJW)
//
//  Numerical version of the Harberger model.  Four household 
//  types differentiated by consumption pattern and income source. 
//----------------------------------------------------------------------//

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

//
//  Parameters describing sectors X and Y
//  
//     Sector X is less capital-intensive
//

decl sx = 0.8 ;
decl dx = 0.4 ; 

decl sy = 0.8 ;
decl dy = 0.6 ;

//
//  Elasticity of substitution in consumption
//

decl sc  = 0.8 ;

//
//  Government preference weight on X 
//

decl dc_g = 0.5 ;

//
//  Factor endowments and preference weight on X
//  by household type.
//
//     Types A and B own capital and labor
//     Types C and D own labor
//
//     Types A and C prefer X
//     Types B and D prefer Y
//

decl lab_a = 200. ;
decl cap_a = 200. ;
decl dc_a  =  0.7 ;

decl lab_b = 200. ;
decl cap_b = 200. ;
decl dc_b  =  0.3 ;

decl lab_c = 200. ;
decl cap_c =   0. ;
decl dc_c  =  0.7 ;

decl lab_d = 200. ;
decl cap_d =   0. ;
decl dc_d  =  0.3 ;

//
//  Taxes and the numeraire price
//

decl tax_kx = 0 ; 
decl tax_ky = 0 ;
decl tax_w  = 0 ;
decl tax_x  = 0 ;
decl tax_y  = 0 ;

decl wage   = 1.0 ;


//
//  These are the endogenous variables:
//

decl rental ; 

decl qx    , qy     ;
decl kx    , ky     ;
decl lx    , ly     ;
decl coc_x , coc_y  ;
decl cost_x, cost_y ; 
decl px    , py     ;

decl pc_a, pc_b, pc_c, pc_d, pc_g ;
decl x_a , x_b , x_c , x_d , x_g  ;
decl y_a , y_b , y_c , y_d , y_g  ;
decl m_a , m_b , m_c , m_d , m_g  ;
decl v_a , v_b , v_c , v_d ;

decl x_agg ;
decl y_agg ;

//
//  These variables are used to store base case values for welfare 
//  calculations and other comparisons.
//

decl tax_kx0, tax_ky0 ;
decl tax_w0 , tax_x0, tax_y0 ;
decl wage0 ;

decl pc_a0, pc_b0, pc_c0, pc_d0, pc_g0 ;
decl m_a0 , m_b0 , m_c0 , m_d0 , m_g0  ;
decl v_a0 , v_b0 , v_c0 , v_d0 ;
decl xbase;

decl rental0 ;
decl px0, py0;
decl qx0, qy0;


//----------------------------------------------------------------------//
//  f(x)
//
//  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(3,1);

   //
   //  Map the guess vector into the model's variables.
   //
   
   rental  = x[0] ;
   qx      = x[1] ;
   qy      = x[2] ;


   //
   //  User cost of capital includes the basic rental price plus the
   //  the sector-specific capital tax
   //
   
   coc_x  = rental + tax_kx ;
   coc_y  = rental + tax_ky ;


   //
   //  Sector X
   //

   cost_x = ( dx*coc_x^(1-sx) + (1-dx)*wage^(1-sx) )^(1/(1-sx)) ;

   kx =   dx   * qx *( cost_x / coc_x )^sx ;
   lx = (1-dx) * qx *( cost_x / wage  )^sx ;

   px = cost_x + tax_x ;


   //
   //  Sector Y
   //

   cost_y = ( dy*coc_y^(1-sy) + (1-dy)*wage^(1-sy) )^(1/(1-sy)) ;
   
   ky =   dy   * qy *( cost_y / coc_y )^sy ;
   ly = (1-dy) * qy *( cost_y / wage  )^sy ;

   py = cost_y + tax_y ;


   //
   //  Household and government income
   //

   m_a = rental*cap_a + wage*(1-tax_w)*lab_a ;
   m_b = rental*cap_b + wage*(1-tax_w)*lab_b ;
   m_c = rental*cap_c + wage*(1-tax_w)*lab_c ;
   m_d = rental*cap_d + wage*(1-tax_w)*lab_d ;

   m_g = tax_kx*kx 
       + tax_ky*ky 
       + tax_w*wage*(lab_a+lab_b+lab_c+lab_d)
       + tax_x*qx 
       + tax_y*qy ;


   //
   //  Consumption price indicies
   //
   
   pc_a = ( dc_a*px^(1-sc) + (1-dc_a)*py^(1-sc) )^(1/(1-sc)) ;
   pc_b = ( dc_b*px^(1-sc) + (1-dc_b)*py^(1-sc) )^(1/(1-sc)) ;
   pc_c = ( dc_c*px^(1-sc) + (1-dc_c)*py^(1-sc) )^(1/(1-sc)) ;
   pc_d = ( dc_d*px^(1-sc) + (1-dc_d)*py^(1-sc) )^(1/(1-sc)) ;

   pc_g = ( dc_g*px^(1-sc) + (1-dc_g)*py^(1-sc) )^(1/(1-sc)) ;


   //
   //  Indirect utility
   //

   v_a = m_a/pc_a ;
   v_b = m_b/pc_b ;
   v_c = m_c/pc_c ;
   v_d = m_d/pc_d ;


   //
   //  Demands for goods
   //

   x_a = (    dc_a  * m_a / pc_a ) * ( pc_a/px )^sc ;
   x_b = (    dc_b  * m_b / pc_b ) * ( pc_b/px )^sc ;
   x_c = (    dc_c  * m_c / pc_c ) * ( pc_c/px )^sc ;
   x_d = (    dc_d  * m_d / pc_d ) * ( pc_d/px )^sc ;
   x_g = (    dc_g  * m_g / pc_g ) * ( pc_g/px )^sc ;

   y_a = ( (1-dc_a) * m_a / pc_a ) * ( pc_a/py )^sc ;
   y_b = ( (1-dc_b) * m_b / pc_b ) * ( pc_b/py )^sc ;
   y_c = ( (1-dc_c) * m_c / pc_c ) * ( pc_c/py )^sc ;
   y_d = ( (1-dc_d) * m_d / pc_d ) * ( pc_d/py )^sc ;
   y_g = ( (1-dc_g) * m_g / pc_g ) * ( pc_g/py )^sc ;

   x_agg = x_a + x_b + x_c + x_d + x_g ;
   y_agg = y_a + y_b + y_c + y_d + y_g ;


   //  Finally, the miss distance is the difference between 
   //  demand and supply in the market for x; omit the labor
   //  market equation because of Walras' Law

   val[0] = kx + ky - (cap_a + cap_b + cap_c + cap_d) ;
   val[1] = qy - y_agg ;
   val[2] = qx - x_agg ;

   return(val);
}


//----------------------------------------------------------------------//
//  printSolution()
//
//  Print out the endogenous variables at the solution.
//----------------------------------------------------------------------//
printSolution()
{
   decl gdp_exp, gdp_inc;
   
   println( " " );
   println( "Walras Law: ", lx+ly-lab_a-lab_b-lab_c-lab_d );
   println( " " );
   println( "tax_kx:  ", tax_kx );
   println( "tax_ky:  ", tax_ky );
   println( "tax_w:   ", tax_w  );
   println( "tax_x:   ", tax_x  );
   println( "tax_y:   ", tax_y  );

   gdp_exp = cost_x*x_agg + cost_y*y_agg ;
   gdp_inc = coc_x*kx + coc_y*ky + wage*(lx+ly) ;
   
   println( " " );
   println( "wage, rental ", wage, rental );
   println( "px, py       ", px, py );
   println( "qx, qy       ", qx, qy );
   println( "gdp exp, inc ", gdp_exp, gdp_inc ); 

   println( " " );
   println( "A:  pc, m, v:  ", pc_a, m_a, v_a );
   println( "B:  pc, m, v:  ", pc_b, m_b, v_b );
   println( "C:  pc, m, v:  ", pc_c, m_c, v_c );
   println( "D:  pc, m, v:  ", pc_d, m_d, v_d );
   println( "Gov pc, m:     ", pc_g, m_g );
}


//----------------------------------------------------------------------//
//  AnalyzePolicy()
//
//  Simulate a policy experiment, print results, calculate and
//  print welfare effects and then reset the policy variables to 
//  their base case values to be ready for the next experiment.
//
//  Also, decompose the EV into price and income parts:
//
//     EV = (Vnew-Vold)*Pold
//     EV = Vnew*Pold - Vold*Pold
//     EV = Vnew*Pold - Vold*Pold + Vnew*Pnew - Vnew*Pnew
//     EV = Vnew*Pold - Mold      + Mnew      - Vnew*Pnew
//     EV = Mnew - Mold - Vnew*(Pnew-Pold)
//     EV = EVm + EVp
//
//     EVm = Mnew - Mold 
//     EVp = - Vnew*(Pnew-Pold)
// 
//----------------------------------------------------------------------//
AnalyzePolicy()
{
   decl ev_a, ev_b, ev_c, ev_d ;
   decl evm_a, evm_b, evm_c, evm_d ;
   decl evp_a, evp_b, evp_c, evp_d ;
   decl totev, excess;
   decl revdiff;

   solveNewton(f,xbase);
   printSolution();

   println( " " );
   println( "Percent Change From Base Case" );
   println( " " );
   println( "wage   ", 100*(wage-wage0)/wage0 );
   println( "rental ", 100*(rental-rental0)/rental0 );
   println( "px     ", 100*(px-px0)/px0 );
   println( "py     ", 100*(py-py0)/py0 );
   println( "qx     ", 100*(qx-qx0)/qx0 );
   println( "qy     ", 100*(qy-qy0)/qy0 );

   ev_a = (v_a - v_a0)*pc_a0;
   ev_b = (v_b - v_b0)*pc_b0;
   ev_c = (v_c - v_c0)*pc_c0;
   ev_d = (v_d - v_d0)*pc_d0;

   evm_a = m_a - m_a0;
   evm_b = m_b - m_b0;
   evm_c = m_c - m_c0;
   evm_d = m_d - m_d0;

   evp_a = -v_a*(pc_a - pc_a0);
   evp_b = -v_b*(pc_b - pc_b0);
   evp_c = -v_c*(pc_c - pc_c0);
   evp_d = -v_d*(pc_d - pc_d0);

   println( " " );
   println( "Welfare Effects, levels" );
   println( " " );
   println( "A:  ev, evm, evp:  ", ev_a, evm_a, evp_a );
   println( "B:  ev, evm, evp:  ", ev_b, evm_b, evp_b );
   println( "C:  ev, evm, evp:  ", ev_c, evm_c, evp_c );
   println( "D:  ev, evm, evp:  ", ev_d, evm_d, evp_d );

   println( " " );
   println( "Welfare Effects, percent of base M" );
   println( " " );
   println( "A:  ev, m, p:  ", 100*ev_a/m_a0, 100*evm_a/m_a0, 100*evp_a/m_a0 );
   println( "B:  ev, m, p:  ", 100*ev_b/m_b0, 100*evm_b/m_b0, 100*evp_b/m_b0 );
   println( "C:  ev, m, p:  ", 100*ev_c/m_c0, 100*evm_c/m_c0, 100*evp_c/m_c0 );
   println( "D:  ev, m, p:  ", 100*ev_d/m_d0, 100*evm_d/m_d0, 100*evp_d/m_d0 );

   totev  = ev_a + ev_b + ev_c + ev_d;
   revdiff= (m_g/pc_g)*pc_g0 - m_g0 ;
   excess = (-totev) - revdiff ;
   
   println( " " );
   println( "Total EV:       ", totev );
   println( "New Revenue:    ", revdiff , " [deflated value]" );
   println( "Excess Burden:  ", excess );
   println( "MEB:            ", excess/revdiff );

   //
   //  Reset policy variables so that the next experiment
   //  will not accidentally carry over some settings from 
   //  this one
   //

   tax_kx = tax_kx0 ; 
   tax_ky = tax_ky0 ;
   tax_w  = tax_w0  ;
   tax_x  = tax_x0  ;
   tax_y  = tax_y0  ;
   wage   = wage0   ;
}


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

main()
{
   decl x;
   decl save, save2 ;

   format( "%9.3f" );

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

   x = ones(3,1);

   println( "" );
   println( "---------------------------" );
   println( " Parameters and Endowments " );
   println( "---------------------------" );
   println( "" );
   println( "X: elast, K weight, L weight: ", sx, dx, 1-dx );
   println( "Y: elast, K weight, L weight: ", sy, dy, 1-dy );
   println( "" );
   println( "Consumption elasticity:       ", sc );
   println( "" );
   println( "A: lab, cap, X weight, Y weight: ", lab_a, cap_a, dc_a, 1-dc_a );
   println( "B: lab, cap, X weight, Y weight: ", lab_b, cap_b, dc_b, 1-dc_b );
   println( "C: lab, cap, X weight, Y weight: ", lab_c, cap_c, dc_c, 1-dc_c );
   println( "D: lab, cap, X weight, Y weight: ", lab_d, cap_d, dc_d, 1-dc_d );
   println( "" );
   println( "Gov X weight, Y weight: ", dc_g, 1-dc_g );

   //
   //  Solve for the base case.
   //

   println( "\f" );
   println( "---------------" );
   println( " (0) Base Case " );
   println( "---------------" );
   println( "" );

   xbase = solveNewton(f,x);
   printSolution();

   //
   //  Save base case variables for welfare calculations
   //  and other comparisions.
   //

   tax_kx0 = tax_kx ; 
   tax_ky0 = tax_ky ;
   tax_w0  = tax_w  ;
   tax_x0  = tax_x  ;
   tax_y0  = tax_y  ;
   wage0   = wage   ;

   pc_a0 = pc_a;
   pc_b0 = pc_b;
   pc_c0 = pc_c;
   pc_d0 = pc_d;
   pc_g0 = pc_g;

   v_a0 = v_a;
   v_b0 = v_b;
   v_c0 = v_c;
   v_d0 = v_d;
 
   m_a0 = m_a;
   m_b0 = m_b;
   m_c0 = m_c;
   m_d0 = m_d;
   m_g0 = m_g;
   
   rental0 = rental ;
   px0     = px ;
   py0     = py ;
   qx0     = qx ;
   qy0     = qy ;

   //
   //  Now run the policy simulations
   //
   
   println( "\f" );
   println( "---------------------" );
   println( " (1) Numeraire Check " );
   println( "---------------------" );
   println( "" );
   
   wage = 2.0*wage ;
   AnalyzePolicy();

   println( "\f" );
   println( "-----------------------" );
   println( " (2) Increase Wage Tax " );
   println( "-----------------------" );
   println( ""  );

   tax_w = tax_w + 0.1 ;
   AnalyzePolicy();

   println( "\f" );
   println( "---------------------------------------" );
   println( " (3) Increase Tax on Capital Used by X " );
   println( "---------------------------------------" );
   println( "" );
   
   tax_kx = tax_kx + 0.55 ;
   AnalyzePolicy();
   
   println( "\f" );
   println( "---------------------------------------" );
   println( " (4) Increase Tax on Capital Used by Y " );
   println( "---------------------------------------" );
   println( "" );
   
   tax_ky = tax_ky + 0.35 ;
   AnalyzePolicy();

   println( "\f" );
   println( "---------------------------------------------------" );
   println( " (5) Equal Increase Tax on Capital Used by X and Y " );
   println( "---------------------------------------------------" );
   println( "" );
   
   tax_kx = tax_kx + 0.2 ;
   tax_ky = tax_ky + 0.2 ;
   AnalyzePolicy();
   
   println( "\f" );
   println( "-----------------------" );
   println( " (6) Increase Tax on X " );
   println( "-----------------------" );
   println( "" );
   
   tax_x = tax_x + 0.145 ;
   AnalyzePolicy();

   println( "\f" );
   println( "-----------------------" );
   println( " (7) Increase Tax on Y " );
   println( "-----------------------" );
   println( "" );
   
   tax_y = tax_y + 0.165 ;
   AnalyzePolicy();

   println( "\f" );
   println( "-------------------------------" );
   println( " (8) Increase Taxes on X and Y " );
   println( "-------------------------------" );
   println( "" );
   
   tax_x = tax_x + 0.072 ;
   tax_y = tax_y + 0.082 ;
   AnalyzePolicy();

}