$TITLE A Sample Implementation of the HHP Decomposition in GAMS

$ontext

Christoph Boehringer
ZEW, Mannheim

Thomas F. Rutherford
University of Colorado

November, 2002

This program illustrates a general purpose procedure for the HHP
decomposition in GAMS models where model solution values (z)
are a implicit function of two or more policy instruments (x).
 
The HHP decomposition attributes changes in the endogenous model
outcome to changes in the exogenous policy variables by evaluating
a line integral from the starting to the ending point.  In this
implementation, all derivatives are calculated numerically.
 
References:

Harrison J., M. Horridge. and K.R.Pearson (2000), Decomposing
Simulation Results with Respect to Exogenous Shocks, Computational
Economics, 15 (3),. 227-249.

Sergey Paltsev, "The Kyoto Protocol: Regional and Sectoral
Contributions to the Carbon Leakage," The Energy Journal, 2001,
22(4), pp. 53-79.

Boehringer, Christoph and Thomas F. Rutherford, "Who should pay how much?
Compensation for International Spillovers from Carbon Abatement
Policies to Developing Countries - A Global CGE Assessment",
Computational Economics, forthcoming.  Currently available at:
../whoshouldpay.pdf

$offtext

*       Retrieve a GAMS library model, GTM, to use for illustration.  

*       We will apply two constraints on network flows and use the 
*       decomposition procedure to determine the contributions of 
*       these two restrictions to demand responses throughout the
*       network.

$call gamslib gtm
$include gtm

*       Declare variables for reporting how restrictions in shipment
*       capacity produce changes in demand:

parameter       decomp          Demand impact associated with network restrictions,
                demand          Demand levels;

*       First consider applying the restrctions one by one and then in concert
*       to illustrate the sensitivity of the decomposition to the order in which
*       the policy shocks are applied.

*       Constrain only us-gulf -> south-west to 2:

x.up("us-gulf","south-west") = 2;
solve gtm maximizing benefit using nlp;
demand(j,"gulf") = d.l(j);

*       Constrain only mid-cont -> south-west to 1:

x.up("us-gulf","south-west") = 3.73;
x.up("mid-cont","south-west") = 1.0;
solve gtm maximizing benefit using nlp;
demand(j,"mid-cont") = d.l(j);

*       Now add a constraint on us-gulf -> south-west:

x.up("us-gulf","south-west") = 2;
solve gtm maximizing benefit using nlp;
demand(j,"both") = d.l(j);

*       Solve the unconstrained case:

x.up("us-gulf","south-west") = 3.73;
x.up("mid-cont","south-west") = 2.3;
solve gtm maximizing benefit using nlp;
demand(j,"ref") = d.l(j);

*       Generate a report for the two alternative decompositions, 
*       based on sequential application of the constraints:

decomp(j,"mid,gulf") = demand(j,"both") - demand(j,"mid-cont");
decomp(j,"gulf,mid") = demand(j,"gulf") - demand(j,"ref");

*       =========================================================================
*       Start the HHP Decomposition Procedure as Implemented in GAMS

*       NB  We use the "_" suffix to avoid naming conflicts with the application
*       and to thereby make code more portable.  There are three places in the
*       procedure which require application-specific statements.

*       --------------------------------------------------------------------
*       1. Application-specific assignments:
*       --------------------------------------------------------------------

*       The following statements relate i_ and j_ to sets in the model.  We
*       need to have these to avoid domain violations:

set i_(i),j_(j);

set     i_              Set of exogenous policy instruments /us-gulf,mid-cont/
        j_              Set of items to be decomposed   / mexico, west-can, ont-quebec, 
                                                  atlantic, new-engl, ny-nj, mid-atl, 
                                                  south-atl, midwest, south-west, central, 
                                                  n-central, west, n-west /;
parameter
        x0_(i_)         Base values of exogenous policy instruments     
                        /us-gulf  3.73, mid-cont  2.3/
        deltax_(i_)     Total change in policy instruments              
                        /us-gulf -1.73, mid-cont -1.3/;


*       --------------------------------------------------------------------
*       Define the precision of the line integral calculation by
*       selecting the number of steps and the numerical differencing
*       interval:

sets    t_              Steps in line integral /t_0*t_100/;

scalar  epsilon_        Differentiation perturbation /0.0001/;


*       Declare other decomposition-specific parameters:

parameters
        dt_                     Step size in integral over t_,
        x_(i_)                  Policy instrument values in the model,
        z_(j_)                  Current value of z,
        z0_(j_,t_)              Values of Z along integral over t_,
        dZdx_(j_,t_,i_)         Values of partial derivatives along integral,
        decomp_(j_,i_)          Decomposition of changes in z,
        handshake_(j_)          Approximation error,
        pctdecomp_(j_,i_)       Decomposition in percentage effects
        isol_                   Solution counter /0/,
        nsol_                   Number of solutions to date;

*   Step size for the line integral:

dt_ = 1 /(card(t_)-1);

*   Initialize policy instrument with base (initial) value before the policy shock 

x_(i_) = x0_(i_);

*       Generate a status report in the title bar if we are operating
*       on Windows NT/2K/XP:

$if %system.filesys%==MSNT file title_ /'title.cmd'/; title_.nd=0; title_.nw=0;

*	This somewhat cryptic command uses the GAMS PUT facility to write
*	a command file which updates the title of the parent window.

*	The CMD file which is written and then executed looks like this:

*	@title Solving case 25 of 303  (8 %% complete)  Start time: 12:28:26, Current time: %time% -- Ctrl-S to pause

*	This command works with Windows NT/2K/XP, but not with Windows 95/98/Me.

$set updatetitle "putclose title_ '@title Solving case ',isol_,' of ',nsol_,'  (',(round(100*isol_/nsol_)),' %% complete)  Start time: %system.time%, Current time: %time% -- Ctrl-S to pause'/; execute 'title.cmd';"

*       Turn off the GAMS column, equation and solution listings:

option limrow=0, limcol=0, solprint=off;

*   Go along the straight line (natural path) for policy instruments in
*   small steps at stepsize dt_

nsol_ = card(t_) * (1+card(i_));
loop(t_,
        isol_ = isol_ + 1;

*       --------------------------------------------------------------------
*       2. Application-specific assignments:
*       --------------------------------------------------------------------

        x.up(i_,"south-west") = x_(i_);

*	Update the title bar with the status prior to the solve:

$if defined title_ %updatetitle%
        solve gtm maximizing benefit using nlp;
        z0_(j_,t_) = d.l(j_);

*   Loop over exogenous policy instruments and compute local derivative

        loop(i_,

                isol_ = isol_ + 1;

*       --------------------------------------------------------------------
*       3. Application-specific assignments:
*       --------------------------------------------------------------------
                x.up(i_,"south-west") = x_(i_) + epsilon_;

*	Update the title bar with the status prior to the solve:

$if defined title_ %updatetitle%
                solve gtm maximizing benefit using nlp;

                x.up(i_,"south-west") = x_(i_);
                z_(j_) = d.l(j_);


*       Numerical differencing evaluates partial derivatives:

                dZdx_(j_,t_,i_) = (z_(j_) - z0_(j_,t_)) / epsilon_;
        );

*       Move input parameter values along the "natural path":

        x_(i_) = x_(i_) + dt_ * deltax_(i_);
);

*   Calculate the absolute change in the model variable Z that is attributable 
*   to the change in the i-th policy instrument

decomp_(j_,i_) = sum(t_, dZdx_(j_,t_,i_) * deltax_(i_) * dt_);

*   Handshake measures the error in adding up
*   (if the error is too large, use more steps to go along the straight line)

handshake_(j_) = sum(i_, decomp_(j_,i_)) 
                - sum(t_$(ord(t_) gt 1), z0_(j_,t_) - z0_(j_,t_-1));

*   Decomposition in percentage effects (share of contribution of the i-th
*   policy instrument in total change of Z)

alias (i_,ii_);

pctdecomp_(j_,i_)$sum(ii_, decomp_(j_,ii_)) 
        = round(100*decomp_(j_,i_)/sum(ii_, decomp_(j_,ii_)), 0); 

display  decomp_, pctdecomp_, handshake_;



*       End of the HHP Decomposition Procedure Implemented in GAMS
*       =========================================================================

*       Generate  a report for comparison with the simple-minded
*       decompositions:

decomp(j_,"HHP") = decomp_(j_,"us-gulf");
decomp(j_,"Tolerance") = handshake_(j_);
display decomp;

$exit


---- 393 PARAMETER decomp_ Decomposition of changes in z us-gulf mid-cont new-engl 0.046 0.017 ny-nj 0.089 0.033 mid-atl 0.066 0.025 south-atl 0.170 0.064 midwest 0.228 0.086 south-west -1.013 -0.295 central -0.089 0.050 n-central -0.126 -0.024 west -0.188 -0.058 n-west -0.045 -0.014 ---- 393 PARAMETER pctdecomp_ Decomposition in percentage effects us-gulf mid-cont new-engl 73.000 27.000 ny-nj 73.000 27.000 mid-atl 73.000 27.000 south-atl 73.000 27.000 midwest 73.000 27.000 south-west 77.000 23.000 central 229.000 -129.000 n-central 84.000 16.000 west 76.000 24.000 n-west 77.000 23.000 ---- 393 PARAMETER handshake_ Approximation error new-engl 0.001, ny-nj 0.002, mid-atl 0.002, south-atl 0.005, midwest 0.006 south-west -0.009, central -0.007, n-central -2.06676E-5, west -0.002, n-west -4.56454E-4
---- 403 PARAMETER decomp Demand impact associated with network restrictions mid,gulf gulf,mid HHP Tolerance new-engl 0.056 0.052 0.046 0.001 ny-nj 0.108 0.100 0.089 0.002 mid-atl 0.081 0.075 0.066 0.002 south-atl 0.209 0.191 0.170 0.005 midwest 0.279 0.167 0.228 0.006 south-west -1.226 -0.937 -1.013 -0.009 central -0.021 -0.139 -0.089 -0.007 n-central -0.139 -0.146 -0.126 -2.06676E-5 west -0.231 -0.171 -0.188 -0.002 n-west -0.055 -0.041 -0.045 -4.56454E-4