### VERSION 4.8.2008 ###

# Housekeeping
reset;
option show_stats 1;


option solver "knitroampl";

model kmax.mod;
data st.dat;
solve;

## OUTPUT ##

printf "\n\n******* Output for kmax ********";
printf "\n K at good steady state: %f\n\n", kmax;

model sttom.mod;

### Generates output matrix of variables evaluated at the Chebyshev nodes ###
set Kspace:= 1..nn;
param ntil {Kspace};
param cytil{Kspace};
param cogtil{Kspace};
param cobtil{Kspace};
param Xtil{Kspace};
param Itil{Kspace};

#printf "param MATRIXtarik : 1    2    3    4  5    6    7    8    9:=\n" >OUTtarik.out;
for {kindex in 1..nn} {
let kinit := K[kindex];
let X:=0.2;
let I:=0.4;
let Rb:=0.1;
let Rg:=0.3;
let cy:=0.2;
let cog:=0.2;
let cob:=0.2;

let nnew:=0.5;
repeat while (abs(nnew-n) > tolerance) 
	{
		if (nnew < 1)
			then {
				let n:=nnew;
				solve;
				let nnew:=(I*(1-gamma)/D)+gamma;
				}
		if (nnew >=1)
			then {
				let n:=1;
				solve;
				let nnew:=I*(1-gamma)/D+gamma;
				if nnew > 1 then let nnew:=1;
				}
	}		

if n==1 then let Rb:=0;
if n==1 then let cob:=0;
if n==1 then let Rg:=alpha*(Q*I)^(alpha-1);
if n==1 then let cy:=(1-alpha)*(kinit^alpha)-I;
if n==1 then let cog:=Rg*Q*I;

printf "%d %f %f %f %f %f %f %f %f %f\n", kindex, K[kindex], n, X, I, Rg, Rb, cog, cob, cy >CWeights.out;
let ntil[kindex]:=n;
let cytil[kindex]:=cy;
let cobtil[kindex]:=cob;
let cogtil[kindex]:=cog;
let Xtil[kindex]:=X;
let Itil[kindex]:=I;
}

####  Generate the Chebyshev weights
param Corder := nn; #actually of order Corder -1, just need

# vectors where the weights will go:
param n_CWeights {k in 1..Corder};
param cy_CWeights {k in 1..Corder};
param cob_CWeights {k in 1..Corder};
param cog_CWeights {k in 1..Corder};
param X_CWeights {k in 1..Corder};
param I_CWeights {k in 1..Corder};

param Tjm{Kspace};

let n_CWeights[1]:=(1/nn)*sum{k in Kspace} ntil[k];
for {mm in 2..Corder} {
    for {j in Kspace}{
    let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial
    }
    let n_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*ntil[x]);
}

let cy_CWeights[1]:=(1/nn)*sum{k in Kspace} cytil[k];
for {mm in 2..Corder} {
    for {j in Kspace}{
    let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial
    }
    let cy_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*cytil[x]);
}

let cob_CWeights[1]:=(1/nn)*sum{k in Kspace} cobtil[k];
for {mm in 2..Corder} {
    for {j in Kspace}{
    let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial
    }
    let cob_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*cobtil[x]);
}

let cog_CWeights[1]:=(1/nn)*sum{k in Kspace} cogtil[k];
for {mm in 2..Corder} {
    for {j in Kspace}{
    let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial
    }
    let cog_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*cogtil[x]);
}

let X_CWeights[1]:=(1/nn)*sum{k in Kspace} Xtil[k];
for {mm in 2..Corder} {
    for {j in Kspace}{
    let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial
    }
    let X_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*Xtil[x]);
}

let I_CWeights[1]:=(1/nn)*sum{k in Kspace} Itil[k];
for {mm in 2..Corder} {
    for {j in Kspace}{
    let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial
    }
    let I_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*Itil[x]);
}

# Display weights of the Chebyshev polynomials 
printf "param :  n_CWeights  cy_CWeights  cog_CWeights cob_CWeights X_CWeights  I_CWeights :=\n" > "CWeights.out";
for {line in 1..nn}{
printf "%d %f %f %f %f %f %f \n", line, n_CWeights[line], cy_CWeights[line],cog_CWeights[line], 
cob_CWeights[line], X_CWeights[line], I_CWeights[line] > "CWeights.out";
}
printf ";" > "CWeights.out";



#####################
# SIMULATIONS
#
#####################

#########################
# Variables definition
#
#########################
var k;
param periods=40; #Number of periods per simulation

param nval;

param outtom;
param meangrowth;
param vargrowth;
param meanoutput;
param varoutput;

param prob;
param captom;
param ytrajectory{1..periods};
param ygrowth{1..(periods-1)};

#additional variables to compute growth rates
var pre ;
fix pre := k;  #intialize
var post;


# Define lower and upper bounds to do simulation
param aux1;
param aux2;
let  aux1 := 0.4*kmax;
let aux2 := 0.6* kmax;

# Number of simulations

param numsim = 20;
param varout{1..numsim};
param meanout{1..numsim};

param kvalues{1..numsim}:= Uniform(aux1,aux2);
display kvalues;

#let k:=kvalues[1];
#display k;



# Call model and data
data CWeights.out;

for {v in 1..numsim}
{
  display v;
  let k := kvalues[v];
  let pre := kvalues[v];
  display k;


   for {ww in 1..periods}{
      # display ww;
      #calculate n for a given k
      let nval:= sum {i in 1..nn} (n_CWeights[i]*cos((i-1)*acos(2*(k-klow)/(khigh-klow)-1)));

      let prob := Uniform(0,1);
      display prob;
      #param prob{nval} = Uniform(0,1);

      # compute capital tomorrow (note that interest is multiplied)
      if prob<=nval then let captom := ( alpha*(  sum {i in 1..nn} (q*X_CWeights[i]*cos((i-1)*acos(2*(k-klow)/(khigh-klow)-1) ) 
				      +Q*I_CWeights[i]*cos(	(i-1)*acos(2*(k-klow)/(khigh-klow)-1))))^(alpha) );
   	       else let captom := ( alpha*(   sum {i in 1..nn}  (q*X_CWeights[i]*cos((i-1)*acos(2*(k-klow)/(khigh-klow)-1))))^(alpha));

      # Redefine capital
      display k;
      display captom;
      printf "\n assigning captom to k and post\n";
      unfix k;
      fix k := captom;
      display k;
      fix post := captom;
      display post;

      # compute output tomorrow (sum young tomorrow to old tomorrow)
      if  prob<=nval then let outtom := (sum {i in 1..nn}  (cy_CWeights[i]+cog_CWeights[i])*cos((i-1)*acos((2*(k-klow)/(kmax-klow)-1))) );
    	       else let  outtom := (sum {i in 1..nn}  (cy_CWeights[i]+cob_CWeights[i])*cos((i-1)*acos((2*(k-klow)/(kmax-klow)-1))) );

      # store values in arrays
      let ytrajectory[ww]:=outtom;
      printf "%d %f \n", ww, ytrajectory[ww]>ytrajectory.out;
      
      #let ygrwth[ww]:=100*(post-pre)/(pre)};
      
      if ww<nn  then let ygrowth[ww]:=100*(post-pre)/(pre);
      if ww<nn then printf "%d %f \n", ww, ygrowth[ww]>ygrowth.out;




      unfix pre;
      fix pre := post; #store for next period
      unfix post;
       display ww;
      };

  # Compute mean output growth
  let meangrowth := (sum {i in 1..(nn-1)} ygrowth[i])/(nn-1);

  #Compute variance of growth

  let vargrowth := (sum {i in 1..(nn-1)} (ygrowth[i]-meangrowth)^2)/(nn-2);

  # Compute variance of output
  let meanoutput := (sum {i in 1..nn} ytrajectory[i] )/(nn);
  let varoutput := (sum {i in 1..(nn-1)} (ygrowth[i]-meanoutput)^2)/(nn-1);

  display ygrowth;
  display ytrajectory;
  display varoutput,vargrowth,meangrowth;
  display varoutput;
  let varout[v] := varoutput;
  let meanout[v] := meanoutput;
};


########################
# Summary statistics
#
########################


param MEAN = (sum {i in 1..numsim} meanout[i])/(numsim);
param meanvar= (sum {i in 1..(numsim-1)} varout[i])/(numsim-1);
param VAR = (sum {i in 1..numsim} (varout[i]-meanvar)^2)/(numsim-2);

display MEAN, VAR;

printf "%f %f\n", ytrajectory>ytrajectory[psi].out;
printf "%f %f\n", ygrowth>ygrowth[psi].out;