$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:
http://debreu.colorado.edu/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);
* Constrain both mid-cont and 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 unconstraint 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:
set t_ Steps in line integral /t_0*t_100/;
scalar epsilon_ Differentiation perturbation /0.0001/;
* Declare other decomposition-specific parameters:
parameter
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,
nsol_ Number of solutions to date;
isol_ = 0;
* 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;
* 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;
$if defined title_ 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';
* --------------------------------------------------------------------
* 2. Application-specific assignments:
* --------------------------------------------------------------------
x.up(i_,"south-west") = x_(i_);
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;
$if defined title_ 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';
* --------------------------------------------------------------------
* 3. Application-specific assignments:
* --------------------------------------------------------------------
x.up(i_,"south-west") = x_(i_) + epsilon_;
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(t_$(ord(t_) gt 1), z0_(j_,t_) - z0_(j_,t_-1))
= Round(100* (sum(i_, decomp_(j_,i_))- sum(t_$(ord(t_) gt 1), z0_(j_,t_) - z0_(j_,t_-1)))
/sum(t_$(ord(t_) gt 1), z0_(j_,t_) - z0_(j_,t_-1)), 0) ;
* 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_,"%_Tol") = handshake_(j_);
display decomp;