R/equilibrium.R

In nopp: Nash Optimal Party Positions

equilibrium <-
function(start, model, data, tolerance=1e-5, max.iter=100, coal=0, alpha=0,
margin=NULL, fixed=NULL, gamma=0, boot=0, MC=0,self.var="self", prox.var="prox", position=NULL, votes=NULL,quadratic=TRUE, conf.level = 0.95){
    ccc <- match.call()
    continue <- TRUE
    iter <- 0
    tmp <- NULL
    k <- length(unique(data$alt))
    if(missing(start))
    start <- sample(data[,self.var], k)
    tmp <- start
    
    basic <- equilibrium.internal(start=start, model=model, data=data, tolerance=tolerance, max.iter=max.iter, coal=coal, alpha=alpha,margin=margin, fixed=fixed, gamma=gamma,
    self.var=self.var, prox.var=prox.var,quadratic=quadratic)
    bootstrap <- NULL
    montecarlo <- NULL
    
    est <- NULL
    mP <- NULL
    
    new.call1 <- as.list(attr(data,"call")) # needed for mlogit.data
    new.call1[[1]] <- NULL
    backup.data <- eval(new.call1$data)
    
    new.call2 <- model$call
    new.call2[[1]] <- NULL
    
    chid <- rownames(backup.data)
    nvoters <- length(chid)
    
    
    if(boot>0){
        cat("\nBootstrapping \n\n")
        n.sim <- boot
        est <- matrix(,n.sim, k)
        mP <- matrix(,n.sim, k)
        pb <- txtProgressBar(min=1, max=n.sim, style=3)
        for(sim in 1:n.sim){
            
            newchid <- sample(chid, nvoters, replace=TRUE)
            #            iidd <- match(data$chid, newchid)
            # iidd <- which(!is.na(iidd))
            # new.election <- data[iidd,]
            tmp.data <- backup.data[newchid,]
            rownames(tmp.data) <- NULL
            mycall <- new.call1
            mycall$data <- as.symbol("tmp.data")
            new.election <- do.call("set.data", mycall)

            mymodelcall <- new.call2
            mymodelcall$data <- as.symbol("new.election")
            mymodelcall <- as.list(mymodelcall)
            new.model <- do.call("mlogit", mymodelcall)

            ans <- try(equilibrium.internal(start=start, model=new.model, data=new.election, tolerance=tolerance, max.iter=max.iter, coal=coal, alpha=alpha, margin=margin, fixed=fixed, gamma=gamma,self.var=self.var, prox.var=prox.var,quadratic=quadratic), "TRUE")
            if(class(ans)[1]=="try-error")
             next()
            #   cat("o")
            est[sim,] <- ans$est
            mP[sim,] <- ans$mP
            setTxtProgressBar(pb, sim)
            #     print(ans)
        }
        close(pb)
        
        est.mean <- apply(est,2,mean)
        est.LB <- apply(est,2,function(x) quantile(x,(1-conf.level)/2))
        est.UB <- apply(est,2,function(x) quantile(x,1-(1-conf.level)/2))
        names(est.mean) <- names(ans$est)
        names(est.LB) <- names(ans$est)
        names(est.UB) <- names(ans$est)
        
        mP.mean <- apply(mP,2,mean)
        mP.LB <- apply(mP,2,function(x) quantile(x,(1-conf.level)/2))
        mP.UB <- apply(mP,2,function(x) quantile(x,1-(1-conf.level)/2))
        names(mP.mean) <- names(ans$mP)
        names(mP.LB) <- names(ans$mP)
        names(mP.UB) <- names(ans$mP)
        
        bootstrap <-  list(est=est, mP=mP, est.mean=est.mean,  est.LB=est.LB, est.UB=est.UB, mP.mean=mP.mean, mP.LB=mP.LB, mP.UB=mP.UB, replications=boot, conf.level = conf.level)
        
        cat("\n")
    }
    
    
    
    if(MC>0){
        cat("\nDoing Monte Carlo \n\n")
        n.sim <- MC
        est <- matrix(,n.sim, k)
        mP <- matrix(,n.sim, k)
        
        IC <- summary(model)$CoefTable
        pb <- txtProgressBar(min=1, max=n.sim, style=3)
        
        for(sim in 1:n.sim){
            
            #newpar <- apply(IC, 1, function(x) rnorm(1, x[1],x[2]))
            newpar <- mvrnorm(1, coef(model), vcov(model))
            
            idx <- match(names(newpar), names(model$coefficients))
            mod1 <- model
            mod1$coefficients[idx] <- newpar
            ans <- try(equilibrium.internal(start=start, model=mod1, data=data, tolerance=tolerance, max.iter=max.iter, coal=coal, alpha=alpha, margin=margin, fixed=fixed, gamma=gamma,
            self.var=self.var, prox.var=prox.var,quadratic=quadratic),"TRUE")
            if(class(ans)[1]=="try-error")
             next()
            #         cat("o")
            est[sim,] <- ans$est
            mP[sim,] <- ans$mP
            setTxtProgressBar(pb, sim)
        }
        close(pb)
        
        est.mean <- apply(est,2,mean)
        est.LB <- apply(est,2,function(x) quantile(x,(1-conf.level)/2))
        est.UB <- apply(est,2,function(x) quantile(x,1-(1-conf.level)/2))
        
        names(est.mean) <- names(ans$est)
        names(est.LB) <- names(ans$est)
        names(est.UB) <- names(ans$est)
        
        mP.mean <- apply(mP,2,mean)
        mP.LB <- apply(mP,2,function(x) quantile(x,(1-conf.level)/2))
        mP.UB <- apply(mP,2,function(x) quantile(x,1-(1-conf.level)/2))
        
        names(mP.mean) <- names(ans$mP)
        names(mP.LB) <- names(ans$mP)
        names(mP.UB) <- names(ans$mP)

        montecarlo <-  list(est=est, mP=mP, est.mean=est.mean, est.LB=est.LB, est.UB=est.UB, mP.mean=mP.mean, mP.LB=mP.LB, mP.UB=mP.UB, replications=MC, conf.level = conf.level)
        cat("\n")
    }
    
    obj <- list(basic=basic, MC=montecarlo,boot=bootstrap, position=position, votes=votes, call=ccc)
    class(obj) <- "nash.eq"
    obj
}


equilibrium.old <-
function(start, model, data, tolerance=1e-5, max.iter=100, coal=0, alpha=0, 
margin=NULL, fixed=NULL, gamma=0, boot=0, MC=0,self.var="self", prox.var="prox", position=NULL, votes=NULL,quadratic=TRUE){
    ccc <- match.call()
    continue <- TRUE
    iter <- 0
    tmp <- NULL
    k <- length(unique(data$alt))
    if(missing(start))
     start <- sample(data[,self.var], k)
    tmp <- start

    basic <- equilibrium.internal(start=start, model=model, data=data, tolerance=tolerance, max.iter=max.iter, coal=coal, alpha=alpha,margin=margin, fixed=fixed, gamma=gamma, 
        self.var=self.var, prox.var=prox.var,quadratic=quadratic)
    bootstrap <- NULL
    montecarlo <- NULL
    
    est <- NULL
    mP <- NULL
    chid <- unique(data$chid)
    nvoters <- length(chid)


    if(boot>0){
        cat("\nBootstrapping \n\n")
        n.sim <- boot
        est <- matrix(,n.sim, k)
        mP <- matrix(,n.sim, k)
        pb <- txtProgressBar(min=1, max=n.sim, style=3)
        for(sim in 1:n.sim){
            
         newchid <- sample(chid, nvoters, replace=TRUE)
         iidd <- match(data$chid, newchid)
         iidd <- which(!is.na(iidd))
         new.election <- data[iidd,]
         ans <- equilibrium.internal(start=start, model=model, data=new.election, tolerance=tolerance, max.iter=max.iter, coal=coal, alpha=alpha, margin=margin, fixed=fixed, gamma=gamma,self.var=self.var, prox.var=prox.var,quadratic=quadratic)
      #   cat("o")
         est[sim,] <- ans$est
         mP[sim,] <- ans$mP  
         setTxtProgressBar(pb, sim)
            #     print(ans)
        }
        close(pb)
        
        est.mean <- apply(est,2,mean)
        est.sd <- apply(est,2,sd)
        names(est.mean) <- names(ans$est)
        names(est.sd) <- names(ans$est)

        mP.mean <- apply(mP,2,mean)
        mP.sd <- apply(mP,2,sd)
        names(mP.mean) <- names(ans$mP)
        names(mP.sd) <- names(ans$mP)
   
        bootstrap <-  list(est.mean=est.mean, est.sd=est.sd, mP.mean=mP.mean, mP.sd=mP.sd,replications=boot)  
        
        cat("\n")
    }  
     
   

    if(MC>0){
        cat("\nDoing Monte Carlo \n\n")
        n.sim <- MC
        est <- matrix(,n.sim, k)
        mP <- matrix(,n.sim, k)
      
        IC <- summary(model)$CoefTable
        pb <- txtProgressBar(min=1, max=n.sim, style=3)
      
        for(sim in 1:n.sim){
        
         newpar <- apply(IC, 1, function(x) rnorm(1, x[1],x[2]))
         idx <- match(names(newpar), names(model$coefficients))
         mod1 <- model
         mod1$coefficients[idx] <- newpar
         ans <- equilibrium.internal(start=start, model=mod1, data=data, tolerance=tolerance, max.iter=max.iter, coal=coal, alpha=alpha, margin=margin, fixed=fixed, gamma=gamma,
             self.var=self.var, prox.var=prox.var,quadratic=quadratic)
#         cat("o")
         est[sim,] <- ans$est
         mP[sim,] <- ans$mP  
         setTxtProgressBar(pb, sim) 
        }
         close(pb)
    
    est.mean <- apply(est,2,mean)
    est.sd <- apply(est,2,sd)
    names(est.mean) <- names(ans$est)
    names(est.sd) <- names(ans$est)

    mP.mean <- apply(mP,2,mean)
    mP.sd <- apply(mP,2,sd)
    names(mP.mean) <- names(ans$mP)
    names(mP.sd) <- names(ans$mP)
     montecarlo <-  list(est.mean=est.mean, est.sd=est.sd, mP.mean=mP.mean, mP.sd=mP.sd,replications=MC)  
        cat("\n")
    }
    
    obj <- list(basic=basic, MC=montecarlo,boot=bootstrap, position=position, votes=votes, call=ccc)
    class(obj) <- "nash.eq"
    obj
}


# ottimizzare spostando margin, gamma, alpha ecc fuori da nash.hip e prepararli
# dentro equilibrium.internal

equilibrium.internal <- function(start, model, data, tolerance=1e-5, max.iter=100, coal=0, alpha=0, margin, fixed=NULL, gamma=0,
    self.var, prox.var,quadratic=TRUE){
    continue <- TRUE
    iter <- 0
    tmp <- start
    names(tmp) <- levels(attr(data,"index")$alt)
    while(continue){
        tmp1 <- nash.hip(tmp, model, data, coal=coal, alpha=alpha, margin=margin, fixed=fixed, gamma=gamma,
            self.var=self.var, prox.var=prox.var,quadratic=quadratic)
        iter <- iter + 1
        #    if(iter==2) continue <- FALSE
        if( (sum(abs(tmp-tmp1$est))<tolerance) | (iter>max.iter) )
        continue <- FALSE
        tmp <- tmp1$est  
     #   cat(".") # print(tmp)
    }
    #  cat("\n")
    tmp1
}


# calculate nash equilibrium
# tutti max voti
# start : posizione del partito
# election: data set agginrato per la distanza elettore-partito
nash.hip <- function(start, model, data, coal=0, alpha=0, margin=NULL, fixed=NULL, gamma=0, self.var, prox.var, quadratic=TRUE){
    if(is.null(coal))
     cat("\nno coalitions")
    
    
    # aggiornamento posizione dei partiti
    #cat("\n")
    #    print(start)
    true.order <- attr(attr(data, "reshapeLong")$varying,"times")
    nvoters <- length(unique(data$chid))
    aa <-  levels(attr(data,"index")$alt) 
    # print(aa)
    newpos <- rep(start[aa], nvoters)
    new.election <- data
    if(quadratic){
     new.election[prox.var] <- -(data[self.var]-newpos)^2
    } else {
     new.election[prox.var] <- -abs(data[self.var]-newpos)
    }
    #print(head(new.election))
    P <- predict(model, newdata=new.election)
    K <- NCOL(P)
  #  print(head(P))
  #  print(str(data))
  #  print(str(new.election))
    
    #print(true.order)
    
    P <- P[, true.order]
    mP <- colMeans(P)
    #  print(mP)
    #    print(colnames(P))
    #print(names(mP))
    
    
    npar <- colnames(P)
    
    
    CC <- rep(0, length(npar) )
    names(CC) <- npar
    if(is.list(coal)){
        C1 <- unlist(coal)
        idx <- match( names(C1), names(CC))
        CC[idx] <- C1
    } else {
        for(i in 1:length(CC))
        CC[i] <- coal[1]
    }
    
    AA <- rep(0, length(npar) ) # default value for alpha=0
    names(AA) <- npar
    
    if(is.list(alpha)){
        A1 <- unlist(alpha)
        idx <- match( names(A1), names(AA))
        AA[idx] <- A1    
    } else {
        for(i in 1:length(AA))
        AA[i] <- alpha[1]
    }
    
    AA[ which(CC==0) ] <- 0  # lonely parties
    
    MM <- vector( mode="list", length(npar))
    names(MM) <- npar
    
    if(is.list(margin)){
        idx <- match( names(margin), names(MM))
        MM[idx] <- margin    
    }
    
    
    # qua sotto formula per max. solo voti del partito
    v2 <- matrix(, nvoters, K )  # Pi*(1-Pi)
    v1 <- matrix(, nvoters, K )  # Pi*(1-Pi)*X
    
    #   for(i in 1:K){
    #    v2[,i] <-  P[,i]*(1-P[,i])
    #    v1[,i] <-  P[,i]*(1-P[,i])*(data$auto[K*(1:nvoters)-(K-1)])  ###
    #}
    
    # formula per coalizioni
    
    for(i in 1:K){
        mm <- 0
        cln <- CC[i] # number of coalition the party belongs to
        idx <- which(CC == cln) # index of parties in the same coalition
        cl <- CC[idx]
        
        if(!is.null(MM[[i]])){
            
            mm <- rowSums(matrix(P[,idx],nvoters,length(idx)))   
            idd <- match(MM[[i]], npar)
            tt <- rowSums(matrix(P[,idd],nvoters,length(idd)))   
            mm <- mm*tt
        }
        
        if(length(idx)>1){
            idx <- idx[-which(names(idx)==names(CC)[i])]
            idx <- as.numeric(idx)
            v2[,i] <- P[,i]*(1-P[,i] - AA[i]*rowSums(matrix(P[, idx],nvoters,length(idx)) ) + mm)
            v1[,i] <- (P[,i]*(1-P[,i] - AA[i]*rowSums(matrix(P[, idx],nvoters,length(idx))) + mm))*(data[K*(1:nvoters)-(K-1),self.var])
            
        } else {
            v2[,i] <-  P[,i]*(1-P[,i] + mm)
            v1[,i] <-  P[,i]*(1-P[,i] + mm)*(data[K*(1:nvoters)-(K-1),self.var])  ###
        }
        
    }
    
    
    
    v1m <- colMeans(v1)
    v2m <- colMeans(v2)
    
    est <- v1m/v2m # nuova posizione partiti
    names(est) <- names(mP)
    
    BB <- numeric(length(npar))
    names(BB) <- npar
    
    if(is.list(gamma)){
        idx <- match( names(gamma), names(BB))
        BB[idx] <- unlist(gamma)    
    } else {
        for(i in 1:length(BB))
        BB[i] <- gamma[1]
    }
    
    
    if(!is.null(fixed)){
        for(i in names(fixed)){
            est[i] <- fixed[[i]] *(1-BB[i]) + est[i]*BB[i]                        
        }
    }
    
    return(list(est=est, mP=mP))
}