```\$TITLE A Sample Implementation of the HHP Decomposition in GAMS

\$ontext

Christoph Boehringer
ZEW, Mannheim

Thomas F. Rutherford

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.

\$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.

*	Set up the solution link to minimize load time:

*       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

```