#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#The fingers game Monte Carlo:
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#tHIS will be a simple strategic probit of this game:
#        /\
#   NA  /  \ A
#      SQ  /\
#     NA  /  \ A
#       Capit War
rm(list=ls())

#link functions
logistic.probA<-function(A,B,param){
    1/(1+exp(param*(B-A)))
}

normal.probA<-function(A,B,param){
    pnorm((A-B)/( (pi^2)/6))   # sqrt(1)
}

#IND==0 means the logistic link
#IND==1 means the normal link
Data.Maker<-function(N,Nvars,b.1a,b.1bc,b.1bd,b.2bc,b.2bd,IND){

    X<-1
    for(i in 2:Nvars){
        X<-cbind(X,rnorm(N,0,1))
    }
    
    B.1a<-rep(b.1a[1],N)
    B.1bc<-rep(b.1bc[1],N)
    B.1bd<-rep(b.1bd[1],N)
    B.2bc<-rep(b.2bc[1],N)
    B.2bd<-rep(b.2bd[1],N)
       
    for(i in 2:Nvars){
        B.1a<-cbind(B.1a,rep(b.1a[i],N))
        B.1bc<-cbind(B.1bc,rep(b.1bc[i],N))
        B.1bd<-cbind(B.1bd,rep(b.1bd[i],N))
        B.2bc<-cbind(B.2bc,rep(b.2bc[i],N))
        B.2bd<-cbind(B.2bd,rep(b.2bd[i],N))
 
    }
    
    u.1a<-apply(B.1a*X,1,sum)
    u.1bc<-apply(B.1bc*X,1,sum)
    u.1bd<-apply(B.1bd*X,1,sum)
    u.2bc<-apply(B.2bc*X,1,sum)
    u.2bd<-apply(B.2bd*X,1,sum)

    if(IND==0){
        p.2c<-logistic.probA(u.2bc,u.2bd,1)
        p.1a<-logistic.probA(u.1a,p.2c*u.1bc+(1-p.2c)*u.1bd,1)
    }
    if(IND==1){
        p.2c<-normal.probA(u.2bc,u.2bd,1)
        p.1a<-normal.probA(u.1a,p.2c*u.1bc+(1-p.2c)*u.1bd,1)
    }

    y.2c<-matrix(NA,nrow=N,ncol=1)
    y.1a<-matrix(NA,nrow=N,ncol=1)
    for(i in 1:N){
        y.2c[i]<-sample(c(1,0),1,prob=c(p.2c[i],(1-p.2c[i])))
        y.1a[i]<-sample(c(1,0),1,prob=c(p.1a[i],(1-p.1a[i])))
    }

    cbind(X,y.1a,y.2c)
}




llik<-function(b,x.1a,x.1bc,x.1bd,x.2bc,x.2bd,x1a.nonzero,x1bc.nonzero,x1bd.nonzero,x2bc.nonzero,x2bd.nonzero,Y,IND){
    
    print(b)
    
    #if nonzero, calculate U1A by combining data with coefficients          
    if(x1a.nonzero==1){
        U1A<-as.matrix(x.1a)%*%as.matrix(b[1:ncol(x.1a)])
        nextcol<-ncol(x.1a)+1
    }
    #if zero, then U1A is zero
    if(x1a.nonzero==0){
        U1A<-as.matrix(x.1a)
        nextcol<-1
    }
    if(x1bc.nonzero==1){
        U1BC<-as.matrix(x.1bc)%*%as.matrix(b[nextcol:(nextcol+ncol(x.1bc)-1)])
        nextcol<-nextcol+ncol(x.1bc)
    }
    if(x1bc.nonzero==0){
        U1BC<-as.matrix(x.1bc)
        nextcol<-nextcol
    }
    if(x1bd.nonzero==1){
        U1BD<-as.matrix(x.1bd)%*%as.matrix(b[nextcol:(nextcol+ncol(x.1bd)-1)])
        nextcol<-nextcol+ncol(x.1bd)
    }
    if(x1bd.nonzero==0){
        U1BD<-as.matrix(x.1bd)
        nextcol<-nextcol
    }
    if(x2bc.nonzero==1){
        U2BC<-as.matrix(x.2bc)%*%as.matrix(b[nextcol:(nextcol+ncol(x.2bc)-1)])
        nextcol<-nextcol+ncol(x.2bc)
    }
    if(x2bc.nonzero==0){
        U2BC<-as.matrix(x.2bc)
        nextcol<-nextcol
    }
    if(x2bd.nonzero==1){
        U2BD<-as.matrix(x.2bd)%*%as.matrix(b[nextcol:(nextcol+ncol(x.2bd)-1)])
        nextcol<-nextcol+ncol(x.2bd)
    }
    if(x2bd.nonzero==0){
        U2BD<-as.matrix(x.2bd)
        nextcol<-nextcol
    }
    
    
    
    if(IND==0){
        p.c<-logistic.probA(U2BC,U2BD,1)
        p.d<-1-p.c
        p.a<-logistic.probA(U1A,(p.c*U1BC+p.d*U1BD),1)
        p.b<-1-p.a
     }
    if(IND==1){
            p.c<-normal.probA(U2BC,U2BD,1)
            p.d<-1-p.c
            p.a<-normal.probA(U1A,(p.c*U1BC+p.d*U1BD),1)
            p.b<-1-p.a
     }
        
    #here are the outcome probabilities
    p.A<-p.a
    p.BC<-p.b*p.c
    p.BD<-p.b*p.d
        
    #now we only need the log-likelihood
    llik.<-Y[,1]*log(p.A)+Y[,2]*log(p.BC)+Y[,3]*log(p.BD)
    sum(llik.)
}




Monte.Carlo<-function(N.sim,N,Nvars,b.1a.T,b.1bc.T,b.1bd.T,b.2bc.T,b.2bd.T,b.1a.Spec,b.1bc.Spec,b.1bd.Spec,b.2bc.Spec,b.2bd.Spec,IND){
    
    B<-matrix(NA,nrow=N.sim,ncol=sum(b.1a.Spec+b.1bc.Spec+b.1bd.Spec+b.2bc.Spec+b.2bd.Spec))
    for(t in 1:N.sim){
        Data<-Data.Maker(N,Nvars,b.1a.T,b.1bc.T,b.1bd.T,b.2bc.T,b.2bd.T,IND)

        Y.a<-Data[,Nvars+1]
        Y.bc<-(1-Y.a)*(Data[,Nvars+2])
        Y.bd<-(1-Y.a)*(1-Data[,Nvars+2])

        Y<-cbind(Y.a,Y.bc,Y.bd)
    
        X<-Data[,1:Nvars]
    
        x1a.nonzero<-1-((sum(b.1a.Spec!=0))==0)*1
        x1bc.nonzero<-1-((sum(b.1bc.Spec!=0))==0)*1
        x1bd.nonzero<-1-((sum(b.1bd.Spec!=0))==0)*1
        x2bc.nonzero<-1-((sum(b.2bc.Spec!=0))==0)*1
        x2bd.nonzero<-1-((sum(b.2bd.Spec!=0))==0)*1
    
        if(x1a.nonzero==1){
            x.1a<-as.matrix(X[,b.1a.Spec*seq(1,Nvars)])
        }else{x.1a<-as.matrix(rep(0,N))}
        if(x1bc.nonzero==1){
            x.1bc<-as.matrix(X[,b.1bc.Spec*seq(1,Nvars)])
        }else{x.1bc<-as.matrix(rep(0,N))}
        if(x1bd.nonzero==1){
            x.1bd<-as.matrix(X[,b.1bd.Spec*seq(1,Nvars)])
        }else{x.1bd<-as.matrix(rep(0,N))}
        if(x2bc.nonzero==1){
            x.2bc<-as.matrix(X[,b.2bc.Spec*seq(1,Nvars)])
        }else{x.2bc<-as.matrix(rep(0,N))}
        if(x2bd.nonzero==1){
            x.2bd<-as.matrix(X[,b.2bd.Spec*seq(1,Nvars)])
        }else{x.2bd<-as.matrix(rep(0,N))}
    
        n.par<-x1a.nonzero*ncol(as.matrix(x.1a))+x1bc.nonzero*ncol(as.matrix(x.1bc))+x1bd.nonzero*ncol(as.matrix(x.1bd))+x2bc.nonzero*ncol(as.matrix(x.2bc))+x2bd.nonzero*ncol(as.matrix(x.2bd))
        stval<-runif(n.par)*0.001
        estimates<-optim(stval,llik,method="BFGS",hessian=TRUE,control=list(fnscale=-1,trace=5,maxit=1500,reltol=1e-15),
        x.1a=x.1a,x.1bc=x.1bc,x.1bd=x.1bd,x.2bc=x.2bc,x.2bd=x.2bd,x1a.nonzero=x1a.nonzero,x1bc.nonzero=x1bc.nonzero,x1bd.nonzero=x1bd.nonzero,x2bc.nonzero=x2bc.nonzero,x2bd.nonzero=x2bd.nonzero,Y=Y,IND=IND)
        B[t,]<-estimates$par
    }
    B
}


#
# Define the true coefficient values.
# First entry is constant.  
# All zeroes is zero.
#
#
true.1a  = c(1,2,0,0,0,0)
true.1bc = c(3,4,0,0,0,0)
true.1bd = c(0,0,0,0,0,0)
true.2bc = c(1,0,0,0,0,2)
true.2bd = c(0,0,0,0,0,0)


#true.1a  = c(1,0,0,0,2,0)
#true.1bc = c(3,4,0,0,0,0)
#true.1bd = c(0,5,0,0,0,0)
#true.2bc = c(1,0,0,0,0,2)
#true.2bd = c(0,0,0,0,0,0)


#
# Define the specification to estimate.
# Entries are 1 or 0 to include or exclude the covariate.
# First entry is constant term.
# All zeroes is zero.
#

spec.1a  = c(1,1,0,0,0,0)
spec.1bc = c(1,1,0,0,0,0)
spec.1bd = c(0,0,0,0,0,0)
spec.2bc = c(1,0,0,0,0,1)
spec.2bd = c(0,0,0,0,0,0)

#spec.1a  = c(1,0,0,0,1,0)
#spec.1bc = c(1,1,0,0,0,0)
#spec.1bd = c(0,1,0,0,0,0)
#spec.2bc = c(1,0,0,0,0,1)
#spec.2bd = c(0,0,0,0,0,0)


# Note that all vectors in true and spec must be the same length


# for debugging

#b.1a.T    = test.1a
#b.1bc.T   = test.1bc
#b.1bd.T   = test.1bc
#b.2bc.T   = test.2bc
#b.2bd.T   = test.2bd
#b.1a.Spec = spec.1a
#b.1bc.Spec = spec.1bc
#b.1bd.Spec = spec.1bd
#b.2bc.Spec = spec.2bc
#b.2bd.Spec = spec.2bd



# We need MC.N to be above 300 or so to avoid difficulties with optim

MC.N.sim=100
MC.N=500
MC.Nvars=6



MC.Result<-Monte.Carlo(MC.N.sim,MC.N,MC.Nvars,true.1a, true.1bc, true.1bd, true.2bc, true.2bd,
         spec.1a, spec.1bc, spec.1bd, spec.2bc, spec.2bd,0)
 

MC.means<-colMeans(MC.Result)
MC.sd<-sd(MC.Result)

True<-c(true.1a,true.1bc,true.1bd,true.2bc,true.2bd)
Spec<-c(spec.1a,spec.1bc,spec.1bd,spec.2bc,spec.2bd)
Spec.Means<-c(spec.1a,spec.1bc,spec.1bd,spec.2bc,spec.2bd)
Spec.SD<-c(spec.1a,spec.1bc,spec.1bd,spec.2bc,spec.2bd)
counter<-1
for(i in 1:length(Spec)){
    if(Spec.Means[i]!=0){
        Spec.Means[i]<-MC.means[counter]
        Spec.SD[i]<-MC.sd[counter]
        counter<-counter+1
    }
    
}

#Table<-matrix(NA,nrow=sum((True+Spec)!=0),ncol=3)
#
#counter<-1
#for(i in 1:length(Spec.Means)){
#    if(True[i]!=0|Spec.Means[i]!=0){
#        Table[counter,1]<-True[i]
#        Table[counter,2]<-Spec.Means[i]
#        Table[counter,3]<-Spec.SD[i]
#        counter<-counter+1
#    }
#}




Result<-matrix(NA,nrow=20,ncol=length(true.1a))
rownames(Result)<-c("1a  True","1a  Est","1a  SD", " ", "1bc True","1bc Est","1bc SD", " ", "1bd True","1bd Est","1bd SD", " ",
"2bc True","2bc Est","2bc SD", " ", "2bd True","2bd Est","2bd SD", " ")

for(i in 1:length(true.1a)){
    for(j in 1:5){
        transl = i + (j-1)*length(true.1a)
        if(True[transl]!=0|Spec.Means[transl]!=0){
            Result[4*(j-1)+1,i]<-True[transl]
            Result[4*(j-1)+2,i]<-Spec.Means[transl]
            Result[4*(j-1)+3,i]<-Spec.SD[transl]
        }
        
    }
    
}



cat("\n")
cat("Results \n")
print(Result)