In skranz/dyngame: Solving stochastic dynamic games with transfers for public perfect equilibria

This RMarkdown document contains the code that replicates the example in Section 4.3 of Goldlücke and Kranz (2017).

Install packages

If you have not yet installed all required R packages, first install dyngame as explained on its Github site:

https://github.com/skranz/dyngame

Cournot example with stochastic water reserves

knitr::opts_chunk$set(tidy=TRUE)

Game definition

library(dyngame)

library(dyngame)

# Some constants used in the game specification

# type of reserve increment
DET = 1  # deterministic reserve increment
UNIF = 2 # uniformely distributed reserve increment

# type of solution
COLL = 1 # collusive solution (optimal dynamic game equilibrium)
MON = 2  # integrated monopoly solution


# a quicker translation function from states vectors to state indices
xv.to.x = function(m,xvm) {
  return( (xvm[,1])*m$xv.dim[[2]]+xvm[,2]+1)
}

# Cournot duopoly with stochastic water reserves. Example from paper
cournot.reserves.game = function(
           x.cap=20, K=20, x.inc=2, delta=2/3,inc.min=0,inc.max=2,
           sol.type=COLL,inc.type=DET,para=NULL) 
{

  # If parameters are given as a list just copy them all
  # into the local environment
  if (!is.null(para)) {
    copy.into.env(source=para)
  }

  # act.fun specifies the action set for a given state xv
  # val = numerical values of actions
  # lab = labels of actions
  act.fun = function(xv,m=NULL) {
    restore.point("act.fun")
    val = list(0:xv[1],0:xv[2])
        lab = val

    list(val=val,lab=lab, i=c(1,2))
  }

  # States can be grouped into sets of states with same
  # action sets. This will speed up computations.
  # In our game each state has a different action set
  x.group = function(xvm,m=NULL) {
    # return a different group index for each state
    x.group = 1:NROW(xvm)
    x.group
  } 

  # Stage game payoff function
  # Specifies the stage game payoffs as a function
  # of the action profile av and state xv
  # The function must be vectorized, i.e. avm and xvm
  # are matrices of action profiles and states
  g.fun = function(avm,xvm,m=NULL) {
    restore.point("g.fun");

        rownames(avm)=rownames(xvm)=NULL

    # compute outputs and prices
        q = avm
    P = K-q[,1]-q[,2]
    P[P<0]=0
    # return profits
    pi = P*q 
    pi
  }

  # Some additional statistics of the solution that we might be interested in.
  # Here we want to know about price, total output, joint profits,
  # consumer surplus, and total welfare 
  extra.sol.fun = function(avm,xvm,m=NULL) {
    q   = avm
    P   = K-q[,1]-q[,2]
    P[P<0]=0
    pi = P*q 
    Q  = q[,1]+q[,2]
    CS = (K-P)*pmin(Q,K)*(1/2)
    Pi = pi[,1]+pi[,2]
    W  = CS+Pi


    mat = cbind(q,P,Q,CS,Pi,W)
    colnames(mat)[1:2]=c("q1","q2")
    return(mat)
  }

  # We do not consider a multistage game
  x.stage.fun = NULL

  # State transitions
  # For a matrix of action profiles and states specifies
  # the matrix of state transitions
  tau.fun = function(avm,xvm,m=NULL) {
    #restore.point("tau.fun")

        rownames(avm)=rownames(xvm)=NULL
    tau = matrix(0,NROW(avm),m$nx)          

        # Firms reserves are restored by a fixed amount x.inc
        if (inc.type==DET) {

      # matrix of resulting state values
            new.xvm = with.floor.and.ceiling(xvm - avm + x.inc ,floor=0,ceiling=x.cap)

            # Translate state values to state indices
      new.x = xv.to.x(m,new.xvm)

      # Set transition probabilities to 1 for resulting states 
            ind.mat = cbind(1:NROW(xvm),new.x)
            tau[ind.mat]=1

        # Alternatively, firms reserves are restored by a uniformely 
    # distributed integer amount between 0 and x.inc
        } else if (inc.type==UNIF) {

      # loop through all combinations of independently
      # distributed restore amount draws
            for (x.add1 in inc.min:inc.max) {
                for (x.add2 in inc.min:inc.max) {
                    x.add = c(x.add1,x.add2)                
                    xvm1 = with.floor.and.ceiling(xvm[,1] - avm[,1] + x.add1 ,floor=0,ceiling=x.cap)
                    xvm2 = with.floor.and.ceiling(xvm[,2] - avm[,2] + x.add2 ,floor=0,ceiling=x.cap)

                    # Translate state values to state indices
          new.x = xv.to.x(m,cbind(xvm1,xvm2))

          # Add probability of restore draw to transition matrix
                    ind.mat = cbind(1:NROW(xvm),new.x)
                    tau[ind.mat]=tau[ind.mat] +1/((inc.max-inc.min+1)^2)
                }
            }   
        }

    tau
  }

  integrated = sol.type == MON

  # return required information for the dynamic game
  list(n=2,
       delta=delta,
       integrated = integrated,
       xv.val = list(0:x.cap,0:x.cap),
       act.fun=act.fun,
       g.fun=g.fun,
       tau.fun=tau.fun,
       x.stage.fun = x.stage.fun,
       x.group=x.group,
       extra.sol.fun=extra.sol.fun)
}

Initializing and solving game

We now solve the game using the parameterizations in our paper.

NO.LABELS = FALSE

# Turn off restore points to speed up computation
set.storing(FALSE)

# Specify parameters 
para = list(
  delta=2/3,
  K=20,
  x.cap=20,
  x.inc=3,
  inc.min=3,
  inc.max=4,
  inc.type=UNIF,
  sol.type=COLL
)

# init game
my.game = cournot.reserves.game(para=para)
mc = init.game(my.game = my.game)

# solve for collusive solution
mc = solve.game(mc)

# Solution table with one row per state
sol = print.sol(mc)
head(sol)

Analyse solution graphically

The code below performs the graphical analysis

# Analyse solution graphically
par(mar=c(5,5,2,2))

# Equilibrium prices
state.levelplot(mc,z=mc$extra.sol.cur[,"P"],
                arrows=FALSE,cuts=10,xrange=c(0,20),zlim=c(7,20),
                main="Prices under Collusion",
                xlab="Reserves firm 1", ylab="Reserves firm 2",
                reverse.colors=TRUE)

# in grey
state.levelplot(mc,z=mc$extra.sol.cur[,"P"],
                arrows=FALSE,cuts=10,xrange=c(0,20),zlim=c(8,20),
                main="Prices under Collusion",
                xlab="Reserves firm 1", ylab="Reserves firm 2",
                reverse.colors=TRUE,col.scheme = "grey")

# Sum of punishment payoffs
V = mc$sol.mat[,"v1"]+mc$sol.mat[,"v2"]
state.levelplot(mc,z=V,arrows=FALSE,cuts=10,xrange=c(0,20),main="Sum of punishment payoffs", xlab="Reserves firm 1", ylab="Reserves firm 2", reverse.colors=!TRUE,col.scheme = "grey")

# Punishment payoffs of player 1      
state.levelplot(mc,z=mc$sol.mat[,"v1"],arrows=FALSE,cuts=10,xrange=c(0,20),xlab="Oil reserves firm 1", ylab="Oil reserves firm 2",reverse.colors=!TRUE,col.scheme = "grey")

skranz/dyngame documentation built on March 27, 2021, 6:03 a.m.