/* Script for solving, simulating, and estimating.
*
*/
reset;
model ag.mod;
data ag.dat;
option solver "./snopt"; # snopt,minos, donlp2, loqo
option snopt_options 'timing 1 outlev=3 iterations=100000';
option randseed 0; # set random seed based on system time
#option solver minos;
#option minos_options 'timing 1 outlev 3';

# set up state
param temp;
for{m in markets, s in 1..nState[m]}
{
 let temp:= (s-1) div nSize;
 for{i in 1..nFirm[m]}
 {
  let state[m,s,i] := temp mod 2;
  let temp := temp div 2;
 }
}   
   
# set data - otherwise it will complain
let{m in markets, t in time} oState[t,m] := 1;
let{m in markets, t in time, i in firms} oAction[t,m,i] := 0;

# initial guesses - equal prob
let{m in markets, s in 1..nState[m], r in 1..nState[m]-1} FP[m,s,r] := 1/nState[m];
let{m in markets, s in 1..nState[m], i in 1..nFirm[m], a in {1}} 
 P[m,i,s,a] := Uniform(0.4,0.6);
let{m in markets, s in 1..nState[m], i in 1..nFirm[m], a in {1}} 
 P[m,i,s,0] := 1-P[m,i,s,1];

display state,size;

# set parameters & solve for equilibrium
fix fixedCost := 1;
fix entryCost := 1;
fix sigma := 1;
fix thetaR := 1;
fix beta := 0.8;
#fix{m in markets, z in 1..nSize-1, w in 1..nSize} sizeTran[m,z,w] := 1/nSize;
fix sizeTran[1,1,1] := 0.75;
fix sizeTran[1,2,2] := 0.75;
#let{m in markets, w in 1..nSize} sizeTran[m,nSize,w] := 1 - sum

drop likelihood;
solve;
printf "\nEnter any number to continue\n";
read temp <-; 

display P;
display state;
display V,size,pi,ep;

printf "\nEnter any number to continue\n";
read temp <-; 
######################################################################
# simulate data
param slast;
param sizenew;
param csum;
param asim{firms};
for{m in markets} 
{
  let slast := ceil(Uniform(1,nState[m]));
  for{t in time} 
  {  
     let oState[t,m] := slast;
     # simulate actions
     for{i in 1..nFirm[m]} 
     {  
       let temp := Uniform(0,1);
       let oAction[t,m,i] := if temp <= P[m,i,slast,1] then 1 else 0;
     }
          
     # update size
     let temp := Uniform(0,1);
     let csum := 0;
     let sizenew := 0;
     repeat while (csum < temp) {
       let sizenew := sizenew + 1;
       let csum := csum + sizeTran[m,sizenew,((slast-1) mod nSize)+1];
     }
     
     # update state
     let slast := 0;
     for{i in 1..nFirm[m]} 
     { 
       let slast := slast + oAction[t,m,i]*2^(i-1);
     }
     let slast := slast*nSize+sizenew
  } # end loop over time
} # end loop over markets

# check simulations
param pState{m in markets, s in states} := 
 1/T * sum{t in time} (if oState[t,m]==s then 1 else 0);
display pState;
param Psim{m in markets, i in firms, s in states, a in actions} := 
 (sum{t in time} (if (oAction[t,m,i]==a && oState[t,m]==s) then 1 else 0)) 
 / (1.0e-307 + sum{t in time} (if oState[t,m]==s then 1 else 0));
display P,Psim;

printf "\nEnter any number to continue\n";
read temp <-;                    

######################################################################
# estimation
unfix entryCost;
unfix fixedCost;
unfix thetaR;
#unfix sigma;
#unfix beta;
unfix sizeTran;

restore likelihood;
# make initial values wrong
#let beta := 0.9;
let fixedCost := 1;
let entryCost := 1;
let thetaR := 1;
# initial guesses - equal prob
#let{m in markets, s in 1..nState[m], r in 1..nState[m]-1} FP[m,s,r] := 1/nState[m];
#let{m in markets, s in 1..nState[m], i in 1..nFirm[m], a in {1}} 
# P[m,i,s,a] := Uniform(0.4,0.6);
#let{m in markets, s in 1..nState[m], i in 1..nFirm[m], a in {1}} 
# P[m,i,s,0] := 1-P[m,i,s,1];

solve;
display entryCost,fixedCost,thetaR,sigma,beta;
printf "\nEnter any number to continue\n";
read temp <-; 
display sizeTran;
display V,P;