# Bernd Beber # Example code for exploring and solving games numerically ############################################################ # plotting utility functions in one dimension # e.g. define u(x)=x^a and graph u(x) against x for different a # define u(x) u <- function(x, a){x^a} # define a set of values for a a <- seq(length=9, from=0, by=.25) # second, set up a screen with multiple rows and columns of graphs par(mfrow=c(3, 3)) # plot utility functions for(i in a) curve(u(x, i), from=0, to=5, ylim=c(0, 5^i), ylab="utility", xlab="value of x", main=paste("utility function for a=", i, sep="")) ############################################################ # plotting utility in two dimensions # e.g. define u(x) = (x^a + y^a)^(1/a) and graph u over x and y as a surface and with contours # compute values for utility function on a grid # set grid resolution n.grid <- 100 y <- x <- seq(0, 5, length=n.grid) # define utility function and parameter u <- function(x, y) (x^a + y^a)^(1/a) a <- .5 # take outer product for values of u on grid z <- outer(x, y, u) # plot surface using "persp" # theta and phi define viewing angle (and ltheta and lphi the illumination angle), r the distance of the eyepoint from the center of the plotting box, expand is the expansion factor for the z-coordinate (here, it shrinks the plot on z) par(mfrow=c(1, 1)) persp(x, y, z, theta=30, phi=30, r=3, expand=.75, border="lightgray", col="lightblue", shade=.7, ltheta=120, lphi=120, ticktype="detailed", zlab="u(x,y)") # plot "contour", for multiple values of a par(mfrow=c(3, 3)) a <- seq(length=9, from=.1, by=.15) for(i in a){ a <- i z <- outer(x, y, u) contour(x, y, z, levels=seq(.05, 2, .15), xlab="values for x", ylab="values for y", main=paste("utility function for a=", i, sep=""), xlim=c(0, 1), ylim=c(0, 1)) } # plot filled contour; a related command is "image" par(mfrow=c(1, 1)) a <- .5 z <- outer(x, y, u) filled.contour(x, y, z) # trellis plotting, with multiple values for a in one plot, using a different example function # use "expand.grid" to create a multidimensional grid require(lattice) rm(x, y, a) n.grid <- 50 plot.grid <- expand.grid(x=seq(-5, 5, length=n.grid), y=seq(-5, 5, length=n.grid), a=c(-1, 4)) attach(plot.grid) u <- function(x, y, a) -((x-a)^2 + (y-a)^2) z <- u(x, y, a) trellis.par.set(theme = col.whitebg()) wireframe(z ~ x * y, groups=a, scales=list(arrows=FALSE), screen=list(z=-50, x=-65), col.groups=colors()[c(26,552)]) wireframe(z ~ x * y, groups=a, scales=list(arrows=FALSE), screen=list(z=-50, x=-65), drape=T, colorkey=T) wireframe(z ~ x * y, groups=a, scales=list(arrows=FALSE), screen=list(z=-50, x=-65), shade=T) # other trellis functions include "levelplot" and "contourplot" ############################################################ # plotting Euclidean preferences for two players in two dimensions # e.g. plot ideal points and indifference curves given ui(x,y) = -(x-ai1)^2 - (y-ai2)^2 for i=1,2 # create an example matrix of players' ideal points ideals <- c(-1, -3, 1, 4) dim(ideals) <- c(2, 2) # utility functions u <- function(x, y, i) -((x - ideals[1,i])^2 + (y - ideals[2,i])^2) u1 <- function(x, y) u(x, y, 1) u2 <- function(x, y) u(x, y, 2) # compute utilities over grid n <- 100 y <- x <- seq(-5, 5, length=n) z <- array(NA, c(n, n, 2)) z[, , 1] <- outer(x, y, u1) z[, , 2] <- outer(x, y, u2) # plot indifference curves contour(x, y, z[, , 1], nlevels=15, xlab="values for x", ylab="values for y", col="gray70") contour(x, y, z[, , 2], nlevels=15, add=T, col="gray30") # add ideal points to plot points(ideals[1, 1], ideals[2, 1], col="gray70", pch=19) points(ideals[1, 2], ideals[2, 2], col="gray30", pch=19) text((t(ideals[1, ])), (t(ideals[2, ])), c("ideal point of player 1", "ideal point of player 2"), pos=1, cex=.7) ############################################################ # plot best response functions for a simple game # e.g let S1 = [0, 1]; S2=[0, 1]; u1 = s1 (a+s2) - s1^2; u2 = s2 (a+s1) - s2^2 # compute utilities over grid, as before n.grid <- 1000 s1 <- s2 <- seq(0, 1, length=n.grid) u1 <- function(s2, s1) s1 * (a + s2) - s1^2 u2 <- function(s1, s2) s2 * (a + s1) - s2^2 a <- .5 u1.values <- outer(s2, s1, u1) u2.values <- outer(s1, s2, u2) # identify best responses best.resp1 <- cbind(s1[max.col(u1.values)], s2) best.resp2 <- cbind(s1, s2[max.col(u2.values)]) # plot 1's best response given any s2, and 2's best response given any s1 # smooth output computed on grid plot(lowess(best.resp1), type="l", xlim=c(0, 1), ylim=c(0, 1), xlab="s1", ylab="s2", col="blue") lines(lowess(best.resp2), col="red") legend(.5, .2, c("best response function for 1", "best response function for 2"), lty=1, col=c("blue", "red"), cex=.7) ############################################################ # iterated elimination of strategies given a payoff matrix # function "iter.elim" for iterated elimination of strictly dominated strategies # inputs: 'row.pay' is the payoff matrix for row player, and 'col.pay' is the payoff matrix for column player iter.elim <- function(row.pay, col.pay) { if(!all(dim(row.pay)==dim(col.pay))) stop("dimensions of matrices don't match") dimnames(col.pay) <- dimnames(row.pay) <- list(paste("r", 1:nrow(row.pay), sep=""), paste("c", 1:ncol(row.pay), sep="")) rm.log <- vector() repeat { rm.row <- logical(nrow(row.pay)) for(i in 1:nrow(row.pay)) { for(j in (1:nrow(row.pay))[-i]) { if(all(row.pay[i,]0, and if each play different strategies they get 0 # then take a 100 row vector, x(0), and randomly assign a 0 or 1 to each element in the vector # randomly select one element from x(0), xj(0), then randomly pull three more elements from x(0), not including xj(0) and identify the optimal response of j to the observed strategies of the other three rows that have been randomly selected # replace this best strategy in x(0) and call the new vector x(1), and do the same thing for each transition from t to t+1 # repeat this e.g. T=1000 times. Graph the share of players playing "1" as t moves from 0 to T # then repeat the whole process e.g. 100 times and graph these 100 series against t # let L be strategy 1, let R be strategy 0 n.pop <- 100; n.rounds <- 1000; n.iter <- 10 n.select <- 1; n.match <- 3 L.pay <- 1; R.pay <- a <- seq(0.5, 2, .5) # set up grid and compute best responses sims <- array(NA, c(n.iter, n.rounds, n.pop, length(R.pay))) for(i in 1:length(R.pay)){ for(j in 1:n.iter){ sims[j, 1, , i] <- round(runif(n.pop)) for(k in 2:n.rounds){ v <- sample(1:n.pop, n.select) w <- sample((1:n.pop)[-v], n.match, replace=T) sims[j, k, , i] <- sims[j, k-1, , i] if(sum((sims[j, k, w, i]==1) * L.pay)==sum((sims[j, k, w, i]==0) * R.pay[i])) { # if both strategies are best responses, a strategy is selected randomly sims[j, k, v, i]=round(runif(n.select)) next } if(sum((sims[j, k, w, i]==1) * L.pay)>sum((sims[j, k, w, i]==0) * R.pay[i])) sims[j, k, v, i]=1 else sims[j, k, v, i]=0 } } } # plot output par(mfrow=c(2, 2)) for(i in 1:length(R.pay)){ plot(c(0, n.rounds), c(0, 1), pch=" ", xlab="Number of rounds", ylab="Share of population playing strategy 1", main=paste("Simulation results for a=", R.pay[i], sep="")) for(j in 1:n.iter){ strat <- vector() for(k in 1:n.rounds){ strat <- c(strat, k, sum(sims[j, k, , i]) / n.pop) } strat <- matrix(strat, ncol=2, byrow=T) lines(strat) } } ############################################################ # finding a stochastically stable state # e.g. take a coordination game in which if both players play L they get 100, if both play R they get 1, and if each play different strategies they get 0 # consider a process with 3 players playing this game # assume that in each period one player is chosen at random and that with probability 1 - e this player chooses his best response to the actions of the other two players, with probability e he chooses the other strategy # define wrapper for multiplying matrix by itself sq <- function(M, t=1) { if(t!=round(t)) { t <- round(t) cat(paste("second argument was rounded to", round(t),"\n")) } if(is.matrix(M)==F) stop("first argument must be matrix") if(ncol(M)!=nrow(M)) stop("first argument must be square matrix") if(t<1) stop("second argument must be greater than or equal to 1") if(t==1) return(M) if(t==2) return(M%*%M) if(t>2) { M.t <- M%*%M for(i in 3:t) M.t <- M.t%*%M return(M.t) } } # set parameter values e <- .001 t <- 1000000 # there are four states in this game, one with nobody playing left, one with one, one with two, and one with three # define transition matrix p <- matrix(c(1-e, e, 0, 0, (1-e)/3, e, 2*(1-e)/3, 0, 0, 2*e/3, (2-e)/3, (1-e)/3, 0, 0, e, (1-e)), 4, 4, byrow=T) sq(p,t) par(mfrow=c(1, 1)) plot(rep(1:4, 4), c(matrix(rep(1:4, 4), 4, byrow=T)), cex=10 * c(matrix(p, 4, byrow=T)), pch=19, col="red", xlim=c(0, 5), ylim=c(0, 5), yaxt="n", xaxt="n", ann=F) # graph for different values of e e <- c(.9, .7, .5, .3, .1, .01) par(mfrow=c(3, 2)) for(i in e) { p.e <- matrix(c(1-i, i, 0, 0, (1-i)/3, i, 2*(1-i)/3, 0, 0, 2*i/3, (2-i)/3, (1-i)/3, 0, 0, i, (1-i)), 4, 4, byrow=T) p.t <- sq(p.e, t) plot(rep(1:4, 4), c(matrix(rep(1:4, 4), 4, byrow=T)), cex=10*c(matrix(p.t, 4, byrow=T)), pch=19, col="red", xlim=c(0, 5), ylim=c(0, 5), yaxt="n", xaxt="n", ann=F) abline(v=c(1:4)) abline(h=c(1:4)) } ############################################################ # function to generate random draws given arbitrary univariate probability density function "pdf" # function "rdraw" to draw from distribution computed on grid (slow for large intervals and large number of decimal places) # inputs: # 'pdf' is the pdf (need not be normalized, but should produce positive values) # 'digits' is number of decimal places (4 by default) # 'from' and 'to' describe the interval over which the pdf is computed (must be bounded, default is unit interval) # 'n' is number of draws rdraw <- function(pdf, from=0, to=1, n=1, digits=4) { grid <- seq(round(from, digits), round(to, digits), 1/10^digits) sample(grid, n, replace=T, prob=pdf(grid)) } # e.g. draws from the uniform distribution rdraw(dunif, n=10) mean(rdraw(dunif, n=10000)) # but preferable to use "runif" runif(n=10) mean(runif(n=10000)) # e.g. consider the following function on interval [1, 2] f <- function(x) 1 - sin(x) rdraw(f, from=1, to=2) # normalizing f for a proper pdf on this interval doesn't change the expected output, because normalization is not needed for "sample" fn <- function(x) (1 - sin(x)) / (1 + cos(2) - cos(1)) rdraw(fn, from=1, to=2) # graphs of the pdf curve(f, xlim=c(1, 2)) curve(fn, xlim=c(1, 2)) integrate(fn, 1, 2) # alternatively, if the cdf (or inverse cdf) is known, draw randomly from uniform and plug those values into the inverse cdf # if only cdf is known, use function "cdf.arg" (or similar optimization) to compute inverse cdf numerically cdf.arg <- function(cdf, cdf.val, from=0, to=1, tol=10^-6) { arg <- rep(NA, length(cdf.val)) for(i in 1:length(cdf.val)) { optim.fct <- function(x) abs(cdf.val[i] - cdf(x)) arg[i] <- optimize(optim.fct, lower=from, upper=to, tol=tol)\$minimum } return(arg) } # e.g. given the same pdf and support as above # define cdf cdf <- function(x) ((x + cos(x)) - (1 + cos(1))) / ((2 + cos(2)) - (1 + cos(1))) # graph unnormalized and normalized cdf curve((x + cos(x)), xlim=c(1, 2)) curve(cdf, xlim=c(1, 2)) # for random draws given pdf, draw randomly from uniform and plug into inverse cdf z <- runif(10000) draws <- cdf.arg(cdf, z, 1, 2) # both methods produce same expected outputs mean(draws) mean(rdraw(f, from=1, to=2, n=10000)) # function "rdraw2" for random draws given arbitrary bivariate probability density function (again, not computationally efficient) rdraw2 <- function(pdf, xlim=c(0, 1), ylim=c(0, 1), n=1, digits=2) { x <- seq(round(xlim[1], digits), round(xlim[2], digits), 1/10^digits) y <- seq(round(ylim[1], digits), round(ylim[2], digits), 1/10^digits) grid <- 1:(length(x) * length(y)) draw.sub <- sample(grid, n, replace=T, prob=outer(x, y, pdf)) draw <- expand.grid(x, y)[draw.sub, ] dimnames(draw) <- list(1:n, c("Var1", "Var2")) return(draw) } # example g <- function(x, y) 2 * x + y rdraw2(g) # graph pdf par(mfrow=c(1, 1)) n.grid <- 100 y <- x <- seq(0, 1, length=n.grid) z <- outer(x, y, g) # plot surface persp(x, y, z, theta=30, phi=30, r=3, expand=.75, border="lightgray", col="lightblue", shade=.7, ltheta=120, lphi=120, ticktype="detailed", zlab="g(x, y)") ############################################################ # finding a counterexample numerically # e.g. find counterexample to the claim that the maximizers of the product are as good at maximizing the sum as the maximizers of the sum themselves, for the following functions # define functions for product and sum of utilities u.prod <- function(x, a) x^a * (1-x)^a u.sum <- function(x, a) x^a + (1-x)^a # set up a grid of values for a and x x.grid <- a.grid <- 10 a <- seq(0, 2, length=a.grid) x <- seq(0, 1, length=x.grid) # compute values of the functions defined above prod.matrix <- outer(x, a, u.prod) sum.matrix <- outer(x, a, u.sum) # now check if claim holds test <- u.sum(x[max.col(t(prod.matrix))],a)==sum.matrix[cbind(max.col(t(sum.matrix)), 1:a.grid)] # plot results plot(a, test, pch=19, cex=2, yaxt="n", ylab="Does the claim hold?", xlab="values of a", main="Finding counterexamples") axis(2, c(0,1), c("No", "Yes")) # using "optimize" instead # set up a vector of values for a a.grid <- seq(0, 2, length=10) # set up a vector to save results of search for counterexamples claim.test <- rep(NA, length(a.grid)) # run a loop to maximize each function and check whether claim holds for (i in 1:length(a.grid)) { u.prod <- function(x, a=a.grid[i]) x^a * (1-x)^a u.sum <- function(x,a=a.grid[i]) x^a + (1-x)^a max.prod <- optimize(u.prod, c(0, 1), maximum=T) max.sum <- optimize(u.sum, c(0, 1), maximum=T) if(u.sum(max.prod\$maximum)==max.sum\$objective) claim.test[i]=T else claim.test[i]=F } # plot results against values of a plot(a.grid, claim.test, pch=19, cex=2, yaxt="n", ylab="Does the claim hold?", xlab="values of a", main="Finding counterexamples") axis(2, c(0, 1), c("No", "Yes"))