$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.
$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:
gtp.solvelink=2;
* 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