/*************** program to solve the model in CCGMV **************/

/**** specification of the output file ****/

new;
format /m1 /ros 10,6;
save path = /programs;
load path = /programs/source;
output file = /programs/result.out reset;


/***** initial conditions (from Campbell and Viceira (1999)) ******/

Msigma = 1/.75 ~ 1 ~ 1/2 ~ 1/4 ~ 1/10 ~ 1/20;
pp = cols(Msigma);

Mgamma = {1 2 4 10 20};
pp2 = cols(Mgamma);

nparam = 4;
ident = eye(nparam);


/* Load files containing parameters from CV */

load t2[9,9] = a0.txt;
t2 = t2[3 5:8,2:3 5:8];

load t3[9,9] = a1.txt;
t3 = t3[3 5:8,2:3 5:8];

load t4[9,9] = b0.txt;
t4 = t4[3 5:8,2:3 5:8];

load t5[9,9] = b1.txt;
t5 = t5[3 5:8,2:3 5:8];

load t6[9,9] = b2.txt;
t6 = t6[3 5:8,2:3 5:8];


/***** parameters *****/

phi = 0.97001;
mu = 0.0687279/4;
Rf = 0.0042278/4;
suc = 0.00225;          /* the conditional sd of shocks to expected return */
su = 0.00925;           /* the unconditional sd of expected return */
rrho = -0.94601;
sv = 0.10769568;        /* the sd of the unexpected market return*/
se = sqrt(1-rrho^2)*sv;
ntheta = rrho*sv/suc;
be = 0.94^0.25;


/***** abscissas and weights for numerical quadrature for n=9 *****/

let ab =    -4.51274586339978    -3.20542900285647     -2.07684797867782
            -1.02325566378913     0.00000000000000      1.02325566378913
             2.07684797867783     3.20542900285647      4.51274586339977;

let wt =    0.00002234584401     0.00278914132123      0.04991640676522
            0.24409750289494     0.40634920634921      0.24409750289494
            0.04991640676522     0.00278914132123      0.00002234584401;

nqp = rows(ab);
abx = ab*suc;
abr = ab*se;


/******* grid for the state variable and ergodic distribution ********/

stepx = 0.0015;
nx = 41;
minx = mu-stepx*(nx-1)/2;
gx0 = seqa(minx,stepx,nx);

/* matrix with the state variables */

const = ones(nx,1);
xx2 = gx0.^2;
xx3 = gx0.^3;
X = const~gx0~xx2~xx3;
X2 = const~gx0~xx2;

/* Compute the ergodic distribution */

ergodic = zeros(1,nx);
ergodic[(nx-1)/2+1] = 1;

epsilon = 1;

do while epsilon>1e-10;
    temperg = zeros(1,nx);
    ind = 1;
    do while ind<=nx;
        if ergodic[ind]>0;
            ind2 = 1;
            do while ind2<=nqp;
                newval = mu+phi*(gx0[ind]-mu)+abx[ind2];
                state = ntoi(newval);
                temperg[state] = temperg[state]+ergodic[ind]*wt[ind2];
                ind2 = ind2+1;
            endo;
        endif;
        ind = ind+1;
    endo;
    temperg = temperg/sumc(temperg');
    epsilon = maxc((abs(ergodic-temperg))');    
    ergodic = temperg;
endo;

save gx0, ergodic;

/* create matrices for saving results */

value = zeros(nx,pp*pp2);
cwratio = zeros(nx,pp*pp2);
portrule = zeros(nparam,pp*pp2);


/***** MAIN LOOP ****/

qq2 = 1;        /* index for gamma */
do while qq2<=pp2;

gam = Mgamma[qq2];

qq1 = 1;        /* index for sigma */
do while qq1<=pp;

    sigma = Msigma[qq1];
    marker = pp*(qq2-1)+qq1;
    print "gamma:" gam "sigma:" sigma;
    

/****** initializing the parameters of the portfolio rule ******/

    param = zeros(nparam,1);
    param[1] = t2[qq2,qq1];    
    param[2] = t3[qq2,qq1];
    print param';

/* Guess the initial value of value function using CV (Eq 24) */
    
    if sigma==1;
        vf0 = const*(1-be);
    else;
        guess = zeros(3,1);
        guess[1] = (t4[qq2,qq1]-sigma*ln(1-be))/(1-sigma); 
        guess[2] = t5[qq2,qq1]/(1-sigma);
        guess[3] = t6[qq2,qq1]/(1-sigma);
        vf0 = exp(X2*guess);
    endif;
    
    ef = ergodic*vf0;
    print ef;


/* SWITCH: 1) GAMMA=SIGMA=1, 2) GAMMA=1,SIGMA!=1, 3) GAMMA!=1 */
    
    if gam==1 and sigma==1;
        vf = vf0;
        cwr = const*(1-be);
    elseif gam==1 and sigma/=1;
        {vf,cwr} = integral(param,vf0);
    else;
    
/***** OPTIMIZING OVER THE COEFFICIENTS ******/

/* initial parameters */

    cc = 1;
    iteration = 1;
    
    if qq2<=2;
        hh = .01;
    elseif qq2<=4;
        hh = .001;
    else;
        hh = .0001;
    endif;

/* first evaluation of the objective function */
    
    if sigma==1;
        fval = intsig(param,vf0);
        xxx1 = const*(1-be);
    else;
        {fval,xxx1} = integral(param,vf0);
    endif;
    
    ef = ergodic*fval;
    print ef;
    
/* intialize variables that hold the best solution */
    
    bestef = ef;
    bestvf = fval;
    bestcwr = xxx1;
    bestparam = param;

/* optimization loop */

    do until (iteration>20 and cc<1e-6) or (iteration>50);
        gradi = zeros(nparam,1);
        hessi = zeros(nparam,nparam);

/* computing the gradient */

        num = 1;

        do while num<=nparam;
            aux1 = param+ident[.,num]*hh;
            aux2 = param-ident[.,num]*hh;

            if sigma==1;
                df1 = intsig(aux1,fval);
                df2 = intsig(aux2,fval);
            else;
                {df1,xxx1} = integral(aux1,fval);
                {df2,xxx1} = integral(aux2,fval);
            endif;
            
            df = (df1-df2)/(2*hh);
            gradi[num] = ergodic*df;
            
            num2 = nparam;
            do while num2>=num;

/* evaluating the function at "x+2h" and "x-2h" to compute the hessian */

                aux1b = aux1+ident[.,num2]*hh;
                aux2b = aux2-ident[.,num2]*hh;
    
                if sigma==1;
                    df1b = intsig(aux1b,fval);
                    df2b = intsig(aux2b,fval);   
                else;
                    {df1b,xxx1} = integral(aux1b,fval);
                    {df2b,xxx1} = integral(aux2b,fval);
                endif;
                
/* computing the hessian */

                auxc = df1b-2*fval+df2b;
                hessi[num,num2] = (ergodic*auxc)/(4*hh^2) ;
                hessi[num2,num] = hessi[num,num2];
                num2 = num2-1;
            endo;
        
            num = num+1;
        endo;
      
        increment = -inv(hessi)*gradi;
        
/* evaluating the objective function at the new set of parameters */

        step = 1;
        nef = 0;
        hu = 0;

        do until nef>=ef or hu>25;
            step = step/1.5;
            tparam = param+step*increment;
            
            if sigma==1;
                vf = intsig(tparam,fval);
                cwr = const*(1-be);
            else;
                {vf,cwr} = integral(tparam,fval);
            endif;
                        
            nef = ergodic*vf;
            hu = hu+1;
        endo;
        
        cc = abs((nef-ef)/ef);

        fval = vf;        
        ef = nef;
        print ef;;
        
        param = tparam;
        print param';

/* save results if optimal thus far */

        if ef>=bestef;
            bestef = ef;
            bestvf = vf;
            bestcwr = cwr;
            bestparam = param;
        endif;
        
        iteration = iteration+1;
    endo;

    vf = bestvf;
    cwr = bestcwr;
    param = bestparam;

    endif;

/* END OF SWITCH */

    
/* print temporary results */

    print "portfolio parameters:";
    print param';
    print "value function & C/W";
    q2 = vf~cwr;
    print q2;

/* save results */

    value[.,marker] = vf;
    cwratio[.,marker] = cwr;
    portrule[.,marker] = param;
    save value, cwratio, portrule;
    
    qq1 = qq1+1;
endo;

qq2 = qq2+1;
endo;

end;


/***********   PROCEDURES   *************/

/* NTOI */
/* procedure that assigns values to corresponding indexes in a grid */

proc ntoi(value);

local index;

index = round((value-minx)/stepx)+1;
index = maxc((minc(index|nx)|1));

retp(index);

endp;


/* INTEGRAL */
/***** procedure that computes the function to be maximized *******/
/* routine that computes CWR, using the Euler equation, for given alfa and V
    then computes the new V and iterates on this until convergence */
    
proc(2) = integral(beta,vf0);

local isgm,numit,maxnumit,minnumit,diff,vf,cwr,v2,
    expect,i,alfa,vvf1,vvf2,A;

isgm = 1-1/sigma;

numit = 1;
maxnumit = 5000;
minnumit = 20;
diff = 1;
vf = vf0;

do until (diff<1e-6 or numit>maxnumit) and (numit>minnumit);
    expect = zeros(nx,1);

    i=1;
    do while i<=nx;
        alfa = X[i,.]*beta;
        expect[i] = vv(vf,i,alfa);
        i = i+1;
    endo;

    if gam==1;
        vvf1 = exp((isgm)*expect);
        vvf2 = exp((1-sigma)*expect);
    else;
        vvf1 = expect.^(isgm/(1-gam));
        vvf2 = expect.^((1-sigma)/(1-gam));
    endif;
    
    A = (((1-be)/be)^sigma)*vvf2;
    cwr = A./(const+A);
    cwr = minc((cwr~0.99*const)');
    cwr = maxc((cwr~1e-10*const)');
    
    v2 = ((1-be)*cwr.^isgm+be*((const-cwr).^isgm).*vvf1).^(1/isgm);
    v2 = minc((v2~10000*const)');
    v2 = maxc((v2~1e-20*const)');
    
    diff = maxc(abs((v2-vf)./vf));
    vf = v2;
    numit = numit+1;
endo;

retp(vf,cwr);

endp;


/* INTSIG */
/***** procedure that computes the function to be maximized for SIGMA=1 *******/

proc intsig(beta,vf0);

local numit,maxnumit,minnumit,diff,vf,cwr,v2,
    expect,i,alfa;

/* routine that, given alfa and V, computes the new V and 
    iterates on this until convergence */

numit = 1;
maxnumit = 5000;
minnumit = 20;
diff = 1;
vf = vf0;

do until (diff<1e-6 or numit>maxnumit) and (numit>minnumit);
    expect = zeros(nx,1);
    
    i=1;
    do while i<=nx;
        alfa = X[i,.]*beta;
        expect[i] = vv(vf,i,alfa);
        i = i+1;
    endo;

    cwr = const*(1-be);
    
    v2 = exp((1-be)*ln(cwr)+be*ln(const-cwr)+be/(1-gam)*ln(expect));
    v2 = minc((v2~10000*const)');
    v2 = maxc((v2~1e-20*const)');

    diff = maxc(abs((v2-vf)./vf));
    vf = v2;
    numit = numit+1;
endo;

retp(vf);

endp;


/* VV */
/* procedure that computes the expectation of the value function
    conditional on one realization of the state variable */

proc vv(v,i,a);

local aux1,aux2,k,j,v0,v1,v2,eps1,num;

v0 = zeros(nqp,1);
v1 = zeros(nqp,1);

k = 1;
do while k<=nqp;

    j = 1;
    do while j<=nqp;
        aux1 = mu+phi*(gx0[i]-mu)+abx[k];
        num = ntoi(aux1);
        
        eps1 = ntheta*abx[k]+abr[j];
        aux2 = a*(exp(gx0[i]+Rf+eps1)-exp(Rf))+exp(Rf);
        aux2 = maxc((aux2~.1)');
        
        if gam==1;
            v0[j] = wt[j]*wt[k]*ln(v[num]*aux2);
        else;
            v0[j] = wt[j]*wt[k]*(v[num]*aux2)^(1-gam);
        endif;
    
        j = j+1;
    endo;        
    
    v1[k] = sumc(v0);
    k= k+1;
endo;
  
v2 = sumc(v1);
retp(v2);

endp;