$ontext
Bellmann value function iteration approach applied to optimal growth model

see. Howitt etal (2002): Using Polynomial Approximation to Solve SDP Problems.

In this version a logarithmic guess of the value function is used as basis
function.

$offtext

option solprint=off, sysout=off
option limrow=0, limcol=0;
$offlisting

*###############################################################################
*                           DATA / SPECIFICATION
*###############################################################################
scalar
         roh     time preference rate    /0.04/
         beta    discount factor
         alpha   capital elasticity      /0.4/
         phi     technology coefficient  /5/
;
beta = 1/(1+roh);
* Set beta to howitt value
beta = 0.9888;

*###############################################################################
*                           COMPUTE VALUE FUNCTION
*###############################################################################
set
         n       number of terms in value function
         /n1*n2/
         k       set of interpolation nodes
         /k1*k10/
;

scalar
         a       lower bound on state    /0.1/
         b       upper bound on state    /5/

parameter
         xhat(k) Chebychev interpolation nodes
         x(k)    converted to "real" interpolation nodes
         tk(k,n) terms of value function approximation
;

*------------------------ COMPUTE INTEPOLATION AND TERMS -----------------------
* Compute interpolation nodes
x(k)     = a + (ord(k)-1)*(b-a)/(card(k)-1);
*display x;

* Determine terms of bisis function
tk(k,"n1") = 1;
tk(k,"n2") = log(x(k));
*display tk;

*--------------------- INITAL GUESS ON VALUE FUNCTION COEFFICIENTS -------------
parameter
         ac(n)   coefficients on polynomial terms
         xcu     current state in value function
;
ac(n) = 1; ac("n1") = 100;


*--------------------- MAXIMIZATION PROBLEM ------------------------------------
Variable
         VAL             value of objective
;

Positive Variable
         XNEXT           next period state variable
         TV              terms of value function
         CONS            consumption
;

equation
         obj     value function
         lom     law of motion
         def_TV1 definition of first polynomial term
         def_TV2 definition of second polynomial term
;

obj..
         VAL                             =E=     log(phi*((xcu+0.0001)**alpha) - XNEXT)
                                                 + beta*sum(n, ac(n)*TV(n))
;

lom..
         XNEXT                           =E=     phi*xcu**alpha - CONS
;


def_TV1(n)$(ord(n) eq 1)..
         TV(n)                           =E=      1
;

def_TV2(n)$(ord(n) eq 2)..
         TV(n)                           =E=      log(XNEXT)
;

model valuefunction /obj, lom, def_TV1, def_TV2/;

xcu   = x("k1");
CONS.LO                  = 1e-5;
XNEXT.L                  = xcu;
XNEXT.LO                 = a;
XNEXT.UP                 = b;
TV.L("n1")               = 1;
TV.L("n2")               = log(xcu);

*------------------------- REGRESSION MODEL ------------------------------------
parameter
         yc      current y value
;

variable
         sqer    squared error sum
         coef    coefficient estimate
;

equation
         objII   objective for regression model
;

objII..
         sqer    =E= sum(k, sqr(yc(k) - sum(n, coef(n)*tk(k,n))))
;

model regress /objII/;
regress.solprint = 2;

*------------------------------- ITERATION -------------------------------------
valuefunction.solvelink=2;
valuefunction.solprint=2;
set
         iter    iterations
         /it1*it1000/
;

scalar
         error   error evaluation        /1/
         conv    convergence criterium   /1e-7/
         stop    stopping flag           /0/
;

parameter
         iterVAL         value function in iterations
         iterNEXT        next period state
         iterTV          polynom terms in iterations
         anew            new coefficients on policnomial terms
         aiterlog        iteration log report for coeffcients
         iterlog         iteration log report other
;

* Iteration loop
loop(iter$(not stop),

*        Interpolation nodes loop
         loop(k,
                 xcu = x(k);
                 XNEXT.UP = phi*(xcu**alpha) - 0.00001;
*                Solve value function
                 solve valuefunction maximizing val using NLP;
*                Store values
                 iterVAL(iter,k)      = VAL.L;
                 iterTV(iter,k,n)     = TV.L(n);
                 iterNEXT(iter)       = XNEXT.L;
         );

*        Obtain new coefficients by regresseion
         yc(k) = iterVAL(iter,k);
         solve regress using NLP minimizing sqer;
         anew(n) = COEF.L(n);

*        Compute error
         error   = sum(n, sqr(anew(n) - ac(n)));

*        Update coefficients in value function
         ac(n)    = anew(n);

*        Write some iteration report
         aiterlog(iter,n)        = ac(n);
         iterlog(iter,"error")   = error;

*        Evaluate error
         stop$(error le conv)    = 1;
);

display iterlog;


*------------------------- COMPARE TO TRUE SOLUTION ----------------------------

parameter
         compare         comparison of true and approximated value function
         const           constant term in true values function
         slope           slope of the value function
;
const = (1/(1-beta))*(
                         log(phi*(1-alpha*beta))
                         + (alpha*beta*log(phi*alpha*beta))/(1-alpha*beta)
                     );

slope = (alpha/(1-alpha*beta));

compare("n1","true")     = const;
compare("n1","computed") = ac("n1");
compare("n2","true")     = slope;
compare("n2","computed") = ac("n2");
display compare;

*###############################################################################
*                           OPTIMAL PATH
*###############################################################################
set
         t       time set
         /t0*t50/
;

parameter
         cap(t)          capital in period t
         k0              initial capital         /0.6/
         compareII       comaprison of policy path
         capbar          steady state capital stock
         solu            solution
;

cap("T0") = k0;
capbar    = (1/(alpha*beta*phi))**(1/(alpha-1));


*--------------------------- TRUE SOLUTION -------------------------------------
variable
         true    true value function value
;

positive variable
         TX      true next period choice
         CONSII  consumption
;

equation
         objIII   true objective
         lomII   law of motion
;

objIII..
         TRUE    =E=     log(CONSII) + const + slope*log(TX)
;

lomII..
         TX      =E=     phi*xcu**alpha - CONSII
;

CONSII.LO = 1e-5;
TX.LO = a ;
TX.UP = b ;
model trueval /objIII, lomII/;
trueval.solvelink = 2;
trueval.solprint = 2;

cap("T0") = k0;
loop(t,
         xcu   = cap(t);
         TX.UP = phi*xcu**alpha - 1e-5;
         solve trueval maximizing TRUE using NLP;
         cap(t+1) = TX.L;
);
compareII(t,"true") = cap(t);

*---------------------------- APPROXIMATED -------------------------------------
cap("T0") = k0;
loop(t,
         xcu = cap(t);
         XNEXT.UP = phi*(xcu**alpha) - 0.00001;
*        Solve value function
         solve valuefunction maximizing val using NLP;
         solu(t,"cap")   = xcu;
         solu(t,"prod")  = phi*xcu**alpha;
         solu(t,"cons")  = phi*xcu**alpha - XNEXT.L;
         solu(t,"inv")   = XNEXT.L;
         solu(t,"obj")   = VAL.L;
         cap(t+1) = XNEXT.L;
);

compareII(t,"approx") = cap(t);

display compareII, capbar, solu;