$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;