$title	A Highly Simplified Stochastic Programming Climate Model

SET	t	Time periods /1*40/;

alias (t,tp);

*	Define decades until we have a breakthrough:

set	sw			States of world /1*10,never/
	transition(t,sw,sw)	State transitions

* Define the event tree by specifying which states generate transitions.
* In this case sw1 defines the root node of the tree.  
*	
		/ 1.1.2 /

loop((t,sw)$transition(t,sw,sw+1),transition(t+1,sw+1,sw+2) =yes;);
display transition;

execute_unload 'probtree.gdx', t, sw, transition;

execute 'gams probtree gdx=treelogic.gdx';

alias (sw,ssw,sow);

set	eq(t,sw)	SW is active in time period T,
	st(t,sw,sow)	SW transitions to SOW in time period T,
	sm(t,sw,sow)	SW is represented by SOW in time period t;

execute_load 'treelogic.gdx', eq, st, sm;

parameter
	sp(t,sw)	Preceding state of world,
	pi(sw)		Period-by-period probability,
	p(t,sw)		Current period-state probability
	var		Variance in probability of event /0.15/;

*	Create a pointer to the preceding state:

loop(st(t,sow,sw),	sp(t+1,sw) =  ord(sw)-ord(sow););

*	Define the probability of an efficiency improvement using a 
*	normal distribution define over a set of time periods:

pi(sw) = exp(-var * sqr(ord(sw)-3)/2);
pi("never") = 0;

*	Assume a 25% chance that there will never be a breakthrough:

pi(sw) = 0.75 * pi(sw)/sum(sow,pi(sow));
*	Define the domain and labels for plotting:

set	tplot(t)	Time periods to plot /1*6/,
	tlabels(t)	/1 2010, 2 2020, 3 2030, 4 2040, 5 2050, 6 2060/

$setglobal domain tplot
$setglobal labels tlabels
$batinclude plot pi

pi("never") = 1 - sum(sw,pi(sw));

*	The current period-state probability:

p(t,sw) = sum(sm(t,sow,sw),pi(sow));

set	t200(t)		Time periods for the first 200 years /1*20/,
	tlbl(t)		Plotting labels for those time periods 
			/5 2050, 10 2100, 15 2150, 20 2200/;

execute_unload 'treeplot.gdx', sw, t200=t, pi, sm, eq,  tlbl=tl;
execute 'gams treeplot --piscale=yes';

parameter
        m0              CO2-equiv concent. 1965 billion tons carbon     /677/,
        tl0             Lower stratum temperature (C) 1965              /.10/,
        t0              Atmospheric temperature (C) 1965                /.2/,
        atret           Marginal atmospheric retention rate             /.64/,
        c1              Coefficient for upper level                     /.226/,
        lam             Climate feedback factor                         /1.41/,
        c3              Coefficient trans upper to lower stratum        /.440/,
        c4              Coeff of transfer for lower level               /.02/,

        r               Rate of social time preference per year         /.03/,
        gl0             Growth rate of population per decade            /.223/,
        dlab            Decline rate of population growth per dec       /.195/,
        deltam          Removal rate carbon per decade                  /.0833/,
        ga0             Initial growth rate for technology per decade   /.15/,
        dela            Decline rate of technology per decade           /.11/,
        gsigma          Growth of sigma per decade                      /-.1168/,
        sig0            CO2-equiv-GWP ratio                             /.519/,
        sigma(t)        Emissions-output ratio,
        L0              1965 world population millions                  /3369/,
        k0              1965 value capital billions 1989 US dollars     /16.03/,
        gamma           Capital elasticity in output                    /.25/,

        a0              Initial level of total factor productivity      /.00963/,
        L(t)            Level of population and labor,
        al(t)           Level of total factor productivity (TFP),
        ga(t)           Growth rate of TFP from 0 to T,
        gl(t)           Growth rate of labor 0 to T,
        gsig(t)         Cumulative improvement of energy efficiency
	ebau(t)		Baseline emissions;

gsig(t) = (gsigma/dela)*(1-EXP(-dela*(ORD(t)-1)));
sigma(t)=sig0*EXP(gsig(t));
gl(t) = (gl0/dlab)*(1-EXP(-dlab*(ORD(t)-1)));
L(t)=L0*EXP(gl(t))*.9;
ga(t)= (ga0/dela)*(1-EXP(-dela*(ORD(t)-1)));
al(t) =a0*EXP(ga(t));

*	Baseline emissions path:

ebau(t) = 10 * sigma(t) * al(t) * (k0*L(t)/L0)**gamma * L(t)**(1-gamma);

parameter
        m(t)            CO2-equiv concentration billion t
        forc(t)         Radiative forcing - W per m2
        forcoth(t)      Exogenous forcings from other greenhouse gases,
        te(t)           Temperature - atmosphere C
        tl(t)           Temperature - lower ocean C
        termv           Terminal value of atmophere
        deltaE          Difference iterval /0.001/;

set	tfirst(t)	The first time period;
tfirst(t) = yes$(ord(t)=1);

*	Non-carbon forcing:

forcoth(t) = 1.42;
m(tfirst) = m0;  
te(tfirst) = t0;  
tl(tfirst) = tl0;

parameter	
	climate		Climate evolution
	eref(t)		Reference emissions
	E(t)		Currently estimated emissions path,
	grad		Temperature gradient,
	teinit		Initial temperature path;

*	Save the baseline emissions path:

*	Write out a "subroutine file" for computing the climate model:

$onecho >climatemodel.gms
loop(t,
	m(t)		= 590 + atret*eref(t) + (1-deltam)*(m(t-1)-590) + m0$tfirst(t);
	forc(t)		= 4.1*(LOG(m(t)/590)/LOG(2)) + forcoth(t);
	te(t)		= te(t-1)+c1*(forc(t-1)-lam*te(t-1)-c3*(te(t-1)-tl(t-1))) + t0$tfirst(t);
	tl(t)		= tl(t-1)+c4*(te(t-1)-tl(t-1)) + tl0$tfirst(t);
);
$offecho

*	Write a second subroutine to compute the gradients:

$onecho >jacobian.gms
	eref(t) = E(t);
$include climatemodel
	teinit(t) = te(t); 
	grad(t,tp) = 0;	
	loop(tp,eref(tp) = eref(tp) + deltaE;
$include climatemodel
	  grad(t,tp)= (te(t)-teinit(t)) / deltaE;
	  te(t) = teinit(t);
	  eref(tp)  = eref(tp) - deltaE;);
$offecho

PARAMETER
	pv(t,sw)	Present value cost
	gradsw		State of world gradients
	tesw		Reference temperature
	erefsw		Reference emissions;

pv(t,sw) = 1/(1+r)**(ord(t)-1);
pv(t,sw)$eq(t-1,sw) = 0.5/(1+r)**(ord(t)-1);
pv(t,"never") = 1/(1+r)**(ord(t)-1);
display pv;


VARIABLES	OBJ		Objective function
		ABATE(t,sw)	Abatement measures;

POSITIVE VARIABLE ABATE;

EQUATIONS	objdef, avetemp;

objdef..  obj =e= sum(eq(t,sw), 
	p(t,sw) * pv(t,sw) * (ABATE(t,sw) 
	+ 10 * sqr(ABATE(t,sw)-ABATE(t-1,sw-sp(t,sw)))/Ebau(t)));

avetemp(t,sw)$eq(t,sw)..	tesw(t,sw) + 
  sum((tp,ssw)$sm(tp,sw,ssw), 
	gradsw(t,tp,sw)*(Ebau(tp)-ABATE(tp,ssw)-erefsw(tp,sw))) =l= 2;

MODEL optabate /objdef, avetemp/;

ABATE.UP(t,sw) = Ebau(t);

PARAMETER	abatelog	Iteration log of Abatement activity;

E(t) = ebau(t);
$include jacobian

tesw(t,sw) = te(t);
gradsw(t,tp,sw) = grad(t,tp);
erefsw(t,sw) = eref(t);
ABATE.L(t,sw) = 0.2 * erefsw(t,sw);

option limrow=1000;

SET	iter	Outer iterations /iter1*iter4/;

loop(iter,

	solve optabate using nlp minimizing obj;
	ABATE.L(t,sw) = sum(sm(t,sw,ssw), ABATE.L(t,ssw));

	abatelog(t,iter) = ABATE.L(t,"never")/ebau(t);

	loop(sw,
	  E(t) = ebau(t) - ABATE.L(t,sw);
$include jacobian
	  tesw(t,sw) = te(t);
	  gradsw(t,tp,sw) = grad(t,tp);
	  erefsw(t,sw) = eref(t);
	);
);

*	Plot the iterative adjustments:

display abatelog;
$libinclude plot abatelog

*	Plot the resulting abatement efforts:

$setglobal gp_opt0 'unset key'
parameter	abateratio(t,sw)	Abatement ratio;
abateratio(t,sw) = ABATE.L(t,sw)/Ebau(t);
$libinclude plot abateratio