R code

Data manipulation:

Appending datasets using text files: (note if you have headers, change accordingly)

#Select files_list

file_names <- list.files(getwd())

file_names <- file_names[grepl(".txt",file_names)]

#Initialize read.csv command

files <- read.csv("sobcov00.txt", header = F, stringsAsFactors = F, sep = "|")

#use lapply to append text files

files <- lapply(file_names, read.csv, header=F, stringsAsFactors = F, sep = "|")

files <- do.call(rbind,files)

files <- setNames(files, c("commyear","statecode", "abb","countycode","countyname","commcode",

"commname", "ip", "ipa","cc", "dt","cl","psc","pepc","pic","uepc", "uic",

"qt","nrq","compacres","la","tpa","sa","ia","lr"))

#write to stata

library(foreign)

write.dta(files, "main_data.dta")

R plots

2. Estimated probability graph: Logit model

#Enter data

dispute.data <- read.table(file="main_data.csv",header=T, sep=",")

head(dispute.data)

#Since multiple specifications were tested: all four model functions are written below, modify accordingly for your own needs:

mod.fit<-glm(formula = dispute ~ totalacres + state, family =

binomial(link = logit), data = dispute.data)

mod.fit2<-glm(formula = dispute ~ totalacres + cdiversity + exp + r + offincome + age_1 + gender +

income.1 + state, family = binomial(link = logit), data = dispute.data)

mod.fit3<-glm(formula = dispute ~ totalacres + cdiversity + exp + totalacres:exp +

r + offincome + age_1 + gender + income.1 + state,

family = binomial(link = logit), data = dispute.data)

mod.fit4<-glm(formula = dispute ~ totalacres + cdiversity + exp + totalacres:exp +

totalacres:r + r + offincome + age_1 + gender + income.1 + state,

family = binomial(link = logit), data = dispute.data)

#Logistic graph

#add first model

linear.pred <- predict(object = mod.fit)

w<- exp(linear.pred)/(1+exp(linear.pred))

head(w)

#add second model prediction

linear.pred2 <- predict(object = mod.fit2)

w2<- exp(linear.pred2)/(1+exp(linear.pred2))

head(w2)

#add third model prediction

linear.pred3 <- predict(object = mod.fit3)

w3 <- exp(linear.pred3)/(1+exp(linear.pred3))

head(w3)

#add third model prediction

linear.pred4 <- predict(object = mod.fit4)

w4 <- exp(linear.pred4)/(1+exp(linear.pred4))

head(w4)

#Create plot with both models

win.graph(width = 7, height = 6, pointsize = 12)

par(bg="white")

plot(x = dispute.data$totalacres, y = w, xlab="Total acres", ylab="Estimated probability",

col = "red",panel.first = grid(col = "gray", lty = "dotted"),

main = "Estimated Likelihood of Commodity Forward Contract Disputes",

ylim=c(0,1), pch=19, lwd=3, axes = F)

axis(side = 1, col = "brown")

axis(side = 2, col = "brown")

points(x = dispute.data$totalacres, y = w2, col = "blue", pch = 24)

points(x = dispute.data$totalacres, y = w3, col = "green", pch = 15)

points(x = dispute.data$totalacres, y = w4, col = "orange", pch = 9)

curve(expr = predict(object = mod.fit,

newdata = data.frame(totalacres=x, state = 1), type = "response"),

col = "red", add = TRUE, lwd = 2,

xlim = c(0, max(dispute.data$totalacres)))

curve(expr = predict(object = mod.fit2,

newdata = data.frame(totalacres=x,cdiversity=1,exp=1,

r=avg_r,offincome=1,age_1=1,gender=1,

income.1=0,state=1), type = "response"),

col = "blue", add = TRUE, lwd = 2, lty = 2,

xlim = c(0, max(dispute.data$totalacres)))

curve(expr = predict(object = mod.fit3,

newdata = data.frame(totalacres= x,cdiversity=1,exp=1,x,

r=avg_r,offincome=1,age_1=1,gender=1,

income.1=0,state=1), type = "response"),

col = "green", add = TRUE, lwd = 2, lty = 3,

xlim = c(0, max(dispute.data$totalacres)))

curve(expr = predict(object = mod.fit4,

newdata = data.frame(totalacres= x,cdiversity=1,exp=1,x,

x,r=avg_r,offincome=1,age_1=1,gender=1,

income.1=0,state=1), type = "response"),

col = "orange", add = TRUE, lwd = 2, lty = 4,

xlim = c(0, max(dispute.data$totalacres)))

legend(x = 2500, y = 0.8, legend = c("Model 1", "Model 2", "Model 3", "Model 4"),

pch = c(12,24,15,9),

col = c("red", "blue", "green", "orange"), bty = "n")

3. Create maps using dplyr and ggplot2

library(data.table)

library(dplyr)

library(ggplot2)

library(maps)

library(mapdata)

library(ggmap)

#Insert Indemnity data for all states

Indemnity_usa <- read.csv("full_indemnity_1.csv")

#Remove potential whitespaces

Indemnity_usa$region <- trimws(Indemnity_usa$region)

Indemnity_usa$subregion <- trimws(Indemnity_usa$subregion)

Indemnity_usa$indemnity <- trimws(Indemnity_usa$indemnity)

#Check structure

str(Indemnity_usa)

#Correct class indemnity

Indemnity_usa$indemnity <- as.numeric(Indemnity_usa$indemnity)

#Create duplicate dataset for adjustment

Indemnity_usa2 <- Indemnity_usa

Indemnity_usa2$indemnity <- Indemnity_usa2$indemnity/1000000

n <- length(Indemnity_usa2$indemnity)

#Adjust max indemnity

for(i in 1:n){

if(Indemnity_usa2$indemnity[i]>=10){

Indemnity_usa2$indemnity[i] <- 10

}

#Check max indemnity

max(Indemnity_usa2$indemnity)

#Create data frame for Indemnity_Map2 and check df

indemnity_map3 <- data.frame(state_names = Indemnity_usa$region, county_names = Indemnity_usa$subregion,

indemnity = Indemnity_usa2$indemnity)

head(indemnity_map3)

#Create county data and check it: to subset further to for eg: a state, use: region == 'ohio'

map.data2 <- data.table(subset(map_data("county")))

head(map.data2)

#Restructure inputs by preference

setkey(map.data2, region, subregion)

#Do the same for indemnity

indemnity_map3 <- data.table(indemnity_map3)

setkey(indemnity_map3, state_names, county_names)

#Create map data for ggplot2 function

map.df3 <- map.data2[indemnity_map3]

str(map.df3)

##Create ggplot2 map by breaking it up

gg <- ggplot(map.df3) + theme_void()

gg <- gg + geom_polygon(aes(x=long, y=lat, group=group, fill=indemnity), color="black")

gg <- gg + coord_map() #Put it in the correct shape

gg <- gg + theme(axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank(),

panel.grid = element_blank(), legend.position=c(0.21,0.1))

gg <- gg + guides(fill = guide_colorbar(barwidth=19.5, title.position = "top", direction="horizontal"))

gg <- gg + scale_fill_gradientn(colors = rev(rainbow(7)),

breaks = c(1,5,9),

"Millions of dollars - red indicates $10 million and more")

gg <- gg + ggtitle("Total Crop Insurance Indemnities - by County - 2017 ")

gg <- gg + theme(plot.title = element_text(margin = margin(t = 20, b = 0 )))

gg <- gg + theme(plot.title = element_text(hjust = 0.5))

gg

Econometric analysis

4. FIML and SUR estimation: this is a challenging code, please be careful when you use this. This is great from manual statistical analysis of econometric models in industrial organization.

#FIML Estimation Using Sem Package:

library(sem)

library(systemfit)

limlRam <- matrix(c(

"consump <- corpProf", "consump_corpProf", NA,

"consump <- corpProfLag", "consump_corpProfLag", NA,

"consump <- wages", "consump_wages", NA,

"invest <- corpProf", "invest_corpProf", NA,

"invest <- corpProfLag", "invest_corpProfLag", NA,

"invest <- capitalLag", "invest_capitalLag", NA,

"privWage <- gnp", "privWage_gnp", NA,

"privWage <- gnpLag", "privWage_gnpLag", NA,

"privWage <- trend", "privWage_trend", NA,

"consump <-> consump", "s11", NA,

"privWage <-> privWage", "s22", NA,

"invest <-> invest", "s33", NA),

ncol = 3, byrow = TRUE)

class(limlRam) <- "mod"

exogVar <- c("corpProf", "corpProfLag", "wages", "capitalLag", "trend", "gnp", "gnpLag")

endogVar <- c("consump", "invest", "privWage")

allVar <- c(exogVar, endogVar)

limlResult <- sem(model = limlRam, S = cov(KleinI[ -1, allVar ]), N = (nrow(KleinI) - 1), fixed.x = exogVar)

attributes(limlResult)

##########################################################################

#M1 Bertrand Nash-Equilibrium

#Q1 = B10 + A11 * P1 + A12 * P2 + B1Y * Y;

#Q2 = B20 + A21 * P1 + A22 * P2 + B2Y * Y;

#P1= -(1/A11) * Q1 + C11 * G + C12 * I1;

#P2= -(1/A22) * Q2 + C22 * I2 ;

FimlRam.M1 <- matrix(c(

"Q1 <-> P1", "A11", -.28,

"Q1 <-> P2", "A12", 0.17,

"Q1 <- Y", "B1Y", NA,

"Q2 <-> P1", "A21", 0.30,

"Q2 <-> P2", "A22", 0.89,

"Q2 <- Y", "B2Y", NA,

"P1 <- G", "C11", NA,

"P1 <- I1", "C12", NA,

"P2 <- I2", "C22", NA,

"Q1 <-> Q1", "s11", NA,

"Q2 <-> Q2", "s22", NA,

"P1 <-> P1", "s33", NA,

"P2 <-> P2", "s44", NA),

ncol = 3, byrow = TRUE)

class(FimlRam.M1) <- "mod.M1"

#Define Endogenous and Exogenous Variables

exogVar <- c("I1","I2","Y","G")

endogVar <- c("Q1","Q2","P1","P2")

allVar <- c(exogVar, endogVar)

#FIML Estimation

FimlResult.M1 <- sem(model = FimlRam.M1, S = cov(data.games[ -1, allVar ]),

N = (nrow(data.games) - 1), fixed.x = exogVar, iterations = 10000)

summary(FimlResult.M1)

#SUR Estimation

#M1 Bertrand Nash-Equilibrium

#Q1 = B10 + A11 * P1 + A12 * P2 + B1Y * Y;

#Q2 = B20 + A21 * P1 + A22 * P2 + B2Y * Y;

#P1= -(1/A11) * Q1 + C11 * G + C12 * I1;

#P2= -(1/A22) * Q2 + C22 * I2 ;

eqQ1 <- Q1~P1+P2+Y

eqQ2 <- Q2~P1+P2+Y

eqP1 <- P1~Q1+G+I1

eqP2 <- P2~Q2+I2

system.M1 <- list(Q1 = eqQ1, Q2 = eqQ2, P1 = eqP1, P2 = eqP2)

surResult.M1 <- systemfit(system.M1, method = "SUR", data = data.games,

methodResidCov = "noDfCor", maxiter = 10000)

summary(surResult.M1)

#Matrices

i <- 1

my.mat <- matrix(c

(1, 2, 3,

4, 5, 6,

7, 8, 9),

byrow=TRUE, ncol=3)

my.mat1 <- matrix(c

(1, 2, 3,

4, 5, 6),

byrow=T, ncol=3)

#Vuong Test for Comparison between model f=Bertrand-Nash and g=Stackelberg in Price US leader:

getwd()

my.Resid <- read.csv(file="Residuals.csv", sep=",", header=T)

my.Resid[1:3,]

attach(my.Resid)

Resid.vcov.M1 <- matrix(c

(141340.7, 112341.7, -106054.2, -15795.06,

112341.7, 149706.8, -171161.4, -9765.94,

-106054.2, -171161.4, 276693.2, -15051.23,

-15795.1, -9765.9, -15051.2, 9199.99),

byrow=T, ncol=4)

Resid.vcov.M2 <- matrix(c

(224404.0, 83985.9, 98425.23, 167855.0,

83985.9, 209602.5, 57092.08, 105156.8,

98425.2, 57092.1, 55582.51, 91786.4,

167855.0, 105156.8, 91786.42, 154177.5),

byrow=T, ncol=4)

Resid.vcov.M3 <- matrix(c

(428950.5, -178577.7, -144785.4, -216183.7,

-178577.7, 84907.8, 36332.7, 71211.6,

-144785.4, 36332.7, 144001.1, 112300.1,

-216183.7, 71211.6, 112300.1, 143522.9),

byrow=T, ncol=4)

Resid.vcov.M4 <- matrix(c

(351095.9, 156217.6, 75758.10, 39683.30,

156217.6, 113919.2, 21270.92, 6333.23,

75758.1, 21270.9, 27026.79, 8883.99,

39683.3, 6333.2, 8883.99, 14556.77),

byrow=T, ncol=4)

Resid.vcov.M5 <- matrix(c

(225095.3, 63347.99, -70112.94, 24244.51,

63348.0, 80565.23, -7798.38, -4657.58,

-70112.9, -7798.38, 30299.99, -4797.93,

24244.5, -4657.58, -4797.93, 12982.94),

byrow=T, ncol=4)

Resid.vcov.M6 <- matrix(c

(1318233, 517302.8, -172070.9, -119290.8,

517303, 223190.3, -63652.0, -52968.2,

-172071, -63652.0, 44510.4, 16151.6,

-119291, -52968.2, 16151.6, 13090.9),

byrow=T, ncol=4)

#M1,M2: f=Bertrand-Nash and g=Stackelberg in Price US leader

df.Resid.M1 <- data.frame(P1= P1[1:72], P2=P2[1:72], Q1=Q1[1:72], Q2=Q2[1:72])

df.Resid.M2 <- data.frame(P1= P1[73:144], P2=P2[73:144], Q1=Q1[73:144], Q2=Q2[73:144])

Resid.Mat.M1 <- as.matrix(df.Resid.M1)

Resid.Mat.M2 <- as.matrix(df.Resid.M2)

Resid.Mat.M1.t <- t(Resid.Mat.M1)

Resid.Mat.M2.t <- t(Resid.Mat.M2)

library(MASS)

Resid.vcov.M1.I <- ginv(Resid.vcov.M1)

Resid.vcov.M2.I <- ginv(Resid.vcov.M2)

rows.Resid.Mat.M1 <- nrow(Resid.Mat.M1)

rows.Resid.Mat.M2 <- nrow(Resid.Mat.M2)

x.M1.M2=matrix(data=NA, nrow=rows.Resid.Mat.M1, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x.M1.M2[j,] <- ((Resid.Mat.M1[j,]%*%Resid.vcov.M1.I%*%Resid.Mat.M1.t[,j])

- (Resid.Mat.M2[j,]%*%Resid.vcov.M2.I%*%Resid.Mat.M2.t[,j]))^2

}

head(x.M1.M2)

normalizer.M1.M2 <- (sqrt(colSums(x.M1.M2)))/2

normalized.LR.M1.M2 <- 2*(-1735.88-(-1765.71))/normalizer.M1.M2

#M1,M3: f=Bertrand-Nash and g=Stackelberg in Price Australia leader

df.Resid.M3 <- data.frame(P1= P1[145:216], P2=P2[145:216], Q1=Q1[145:216], Q2=Q2[145:216])

Resid.Mat.M3 <- as.matrix(df.Resid.M3)

Resid.Mat.M3.t <- t(Resid.Mat.M3)

Resid.vcov.M3.I <- ginv(Resid.vcov.M3)

rows.Resid.Mat.M3 <- nrow(Resid.Mat.M3)

x1.M1.M3=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M1.M3[j,] <- ((Resid.Mat.M1[j,]%*%Resid.vcov.M1.I%*%Resid.Mat.M1.t[,j])

- (Resid.Mat.M3[j,]%*%Resid.vcov.M3.I%*%Resid.Mat.M3.t[,j]))^2

}

head(x1.M1.M3)

normalizer.M1.M3 <- (sqrt(colSums(x1.M1.M3)))/2

normalized.LR.M1.M3 <- 2*(-1735.88-(-1736.15))/normalizer.M1.M3

normalized.LR.M1.M3

#M1,M4: f=Bertrand-Nash and g=Cournot Nash

df.Resid.M4 <- data.frame(P1= P1[217:288], P2=P2[217:288], Q1=Q1[217:288], Q2=Q2[217:288])

Resid.Mat.M4 <- as.matrix(df.Resid.M4)

Resid.Mat.M4.t <- t(Resid.Mat.M4)

Resid.vcov.M4.I <- ginv(Resid.vcov.M4)

rows.Resid.Mat.M4 <- nrow(Resid.Mat.M4)

x1.M1.M4=matrix(data=NA, nrow=rows.Resid.Mat.M4, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M1.M4[j,] <- ((Resid.Mat.M1[j,]%*%Resid.vcov.M1.I%*%Resid.Mat.M1.t[,j])

- (Resid.Mat.M4[j,]%*%Resid.vcov.M4.I%*%Resid.Mat.M4.t[,j]))^2

}

head(x1.M1.M4)

normalizer.M1.M4 <- (sqrt(colSums(x1.M1.M4)))/2

normalized.LR.M1.M4 <- 2*(-1735.88-(-1786.33))/normalizer.M1.M4

normalized.LR.M1.M4

#M1,M5: f=Bertrand-Nash and g=Stackelberg Quantity, US Leader.

df.Resid.M5 <- data.frame(P1= P1[289:360], P2=P2[289:360], Q1=Q1[289:360], Q2=Q2[289:360])

Resid.Mat.M5 <- as.matrix(df.Resid.M5)

Resid.Mat.M5.t <- t(Resid.Mat.M5)

Resid.vcov.M5.I <- ginv(Resid.vcov.M5)

rows.Resid.Mat.M5 <- nrow(Resid.Mat.M5)

x1.M1.M5=matrix(data=NA, nrow=rows.Resid.Mat.M4, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M1.M5[j,] <- ((Resid.Mat.M1[j,]%*%Resid.vcov.M1.I%*%Resid.Mat.M1.t[,j])

- (Resid.Mat.M5[j,]%*%Resid.vcov.M5.I%*%Resid.Mat.M5.t[,j]))^2

}

head(x1.M1.M5)

normalizer.M1.M5 <- (sqrt(colSums(x1.M1.M5)))/2

normalized.LR.M1.M5 <- 2*(-1735.88-(-1765.97))/normalizer.M1.M5

normalized.LR.M1.M5

#M1,M6: f=Bertrand-Nash and g=Stackelberg Quantity, Australia Leader.

df.Resid.M6 <- data.frame(P1= P1[361:432], P2=P2[361:432], Q1=Q1[361:432], Q2=Q2[361:432])

Resid.Mat.M6 <- as.matrix(df.Resid.M6)

Resid.Mat.M6.t <- t(Resid.Mat.M6)

Resid.vcov.M6.I <- ginv(Resid.vcov.M6)

rows.Resid.Mat.M5 <- nrow(Resid.Mat.M6)

x1.M1.M6=matrix(data=NA, nrow=rows.Resid.Mat.M4, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M1.M6[j,] <- ((Resid.Mat.M1[j,]%*%Resid.vcov.M1.I%*%Resid.Mat.M1.t[,j])

- (Resid.Mat.M6[j,]%*%Resid.vcov.M6.I%*%Resid.Mat.M6.t[,j]))^2

}

head(x1.M1.M6)

normalizer.M1.M6 <- (sqrt(colSums(x1.M1.M5)))/2

normalized.LR.M1.M6 <- 2*(-1735.88-(-1733.8))/normalizer.M1.M6

normalized.LR.M1.M6

#M2,M3: f=Stackelberg Price US Leader and g=Stackelberg Price, Oz Leader.

x1.M2.M3=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M2.M3[j,] <- ((Resid.Mat.M2[j,]%*%Resid.vcov.M2.I%*%Resid.Mat.M2.t[,j])

- (Resid.Mat.M3[j,]%*%Resid.vcov.M3.I%*%Resid.Mat.M3.t[,j]))^2

}

head(x1.M2.M3)

normalizer.M2.M3 <- (sqrt(colSums(x1.M2.M3)))/2

normalized.LR.M2.M3 <- 2*(-1765.71-(-1736.15))/normalizer.M2.M3

normalized.LR.M2.M3

#M2,M4: f=Stackelberg Price US Leader and g=Cournot-Nash.

x1.M2.M4=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M2.M4[j,] <- ((Resid.Mat.M2[j,]%*%Resid.vcov.M2.I%*%Resid.Mat.M2.t[,j])

- (Resid.Mat.M4[j,]%*%Resid.vcov.M4.I%*%Resid.Mat.M4.t[,j]))^2

}

head(x1.M2.M4)

normalizer.M2.M4 <- (sqrt(colSums(x1.M2.M4)))/2

normalized.LR.M2.M4 <- 2*(-1765.71-(-1786.33))/normalizer.M2.M4

normalized.LR.M2.M4

#M2,M5: f=Stackelberg Price US Leader and g=Stackelberg Quantity, US Leader.

x1.M2.M5=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M2.M5[j,] <- ((Resid.Mat.M2[j,]%*%Resid.vcov.M2.I%*%Resid.Mat.M2.t[,j])

- (Resid.Mat.M5[j,]%*%Resid.vcov.M5.I%*%Resid.Mat.M5.t[,j]))^2

}

head(x1.M2.M5)

normalizer.M2.M5 <- (sqrt(colSums(x1.M2.M5)))/2

normalized.LR.M2.M5 <- 2*(-1765.71-(-1765.97))/normalizer.M2.M5

normalized.LR.M2.M5

#M2,M6: f=Stackelberg Price US Leader and g=Stackelberg Quantity, Oz Leader.

x1.M2.M6=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M2.M6[j,] <- ((Resid.Mat.M2[j,]%*%Resid.vcov.M2.I%*%Resid.Mat.M2.t[,j])

- (Resid.Mat.M6[j,]%*%Resid.vcov.M6.I%*%Resid.Mat.M6.t[,j]))^2

}

head(x1.M2.M6)

normalizer.M2.M6 <- (sqrt(colSums(x1.M2.M6)))/2

normalized.LR.M2.M6 <- 2*(-1765.71-(-1733.8))/normalizer.M2.M6

normalized.LR.M2.M6

#M3,M4: f=Stackelberg Price Oz Leader and g=Cournot-Nash.

x1.M3.M4=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M3.M4[j,] <- ((Resid.Mat.M3[j,]%*%Resid.vcov.M3.I%*%Resid.Mat.M3.t[,j])

- (Resid.Mat.M4[j,]%*%Resid.vcov.M4.I%*%Resid.Mat.M4.t[,j]))^2

}

head(x1.M3.M4)

normalizer.M3.M4 <- (sqrt(colSums(x1.M3.M4)))/2

normalized.LR.M3.M4 <- 2*(-1736.15-(-1786.33))/normalizer.M3.M4

normalized.LR.M3.M4

#M3,M5: f=Stackelberg Price Oz Leader and g=Stackelberg Quantity US Leader. .

x1.M3.M5=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M3.M5[j,] <- ((Resid.Mat.M3[j,]%*%Resid.vcov.M3.I%*%Resid.Mat.M3.t[,j])

- (Resid.Mat.M5[j,]%*%Resid.vcov.M5.I%*%Resid.Mat.M5.t[,j]))^2

}

head(x1.M3.M5)

normalizer.M3.M5 <- (sqrt(colSums(x1.M3.M5)))/2

normalized.LR.M3.M5 <- 2*(-1736.15-(-1765.97))/normalizer.M3.M5

normalized.LR.M3.M5

#M3,M6: f=Stackelberg Price Oz Leader and g=Stackelberg Quantity Oz Leader. .

x1.M3.M6=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M3.M6[j,] <- ((Resid.Mat.M3[j,]%*%Resid.vcov.M3.I%*%Resid.Mat.M3.t[,j])

- (Resid.Mat.M6[j,]%*%Resid.vcov.M6.I%*%Resid.Mat.M6.t[,j]))^2

}

head(x1.M3.M6)

normalizer.M3.M6 <- (sqrt(colSums(x1.M3.M6)))/2

normalized.LR.M3.M6 <- 2*(-1736.15-(-1733.8))/normalizer.M3.M6

normalized.LR.M3.M6

#M4,M5: f=Cournot-Nash and g=Stackelberg Quantity US Leader. .

x1.M4.M5=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M4.M5[j,] <- ((Resid.Mat.M4[j,]%*%Resid.vcov.M4.I%*%Resid.Mat.M4.t[,j])

- (Resid.Mat.M5[j,]%*%Resid.vcov.M5.I%*%Resid.Mat.M5.t[,j]))^2

}

head(x1.M4.M5)

normalizer.M4.M5 <- (sqrt(colSums(x1.M4.M5)))/2

normalized.LR.M4.M5 <- 2*(-1786.33-(-1765.97))/normalizer.M4.M5

normalized.LR.M4.M5

#M4,M6: f=Cournot-Nash and g=Stackelberg Quantity Oz Leader. .

x1.M4.M6=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M4.M6[j,] <- ((Resid.Mat.M4[j,]%*%Resid.vcov.M4.I%*%Resid.Mat.M4.t[,j])

- (Resid.Mat.M6[j,]%*%Resid.vcov.M6.I%*%Resid.Mat.M6.t[,j]))^2

}

head(x1.M4.M6)

normalizer.M4.M6 <- (sqrt(colSums(x1.M4.M6)))/2

normalized.LR.M4.M6 <- 2*(-1786.33-(-1733.8))/normalizer.M4.M6

normalized.LR.M4.M6

#M5,M6: f=Stackelberg Quantity, US Leader and g=Stackelberg Quantity Oz Leader. .

x1.M5.M6=matrix(data=NA, nrow=rows.Resid.Mat.M3, ncol=1)

for(j in 1:rows.Resid.Mat.M1){

x1.M5.M6[j,] <- ((Resid.Mat.M4[j,]%*%Resid.vcov.M5.I%*%Resid.Mat.M5.t[,j])

- (Resid.Mat.M6[j,]%*%Resid.vcov.M6.I%*%Resid.Mat.M6.t[,j]))^2

}

head(x1.M5.M6)

normalizer.M5.M6 <- (sqrt(colSums(x1.M5.M6)))/2

normalized.LR.M5.M6 <- 2*(-1765.97-(-1733.8))/normalizer.M5.M6

normalized.LR.M5.M6

#Derive VCOV: M1 and M2

df.Resid.M1.new <- data.frame(Q1=Q1[1:72], Q2=Q2[1:72], P1= P1[1:72], P2=P2[1:72])

df.Resid.M2.new <- data.frame(Q1=Q1[73:144], Q2=Q2[73:144], P1= P1[73:144], P2=P2[73:144])

Resid.Mat.M1 <- as.matrix(df.Resid.M1.new)

Resid.Mat.M2 <- as.matrix(df.Resid.M2.new)

Resid.Mat.M1.t <- t(Resid.Mat.M1)

Resid.Mat.M2.t <- t(Resid.Mat.M2)

new.vcov.M1 <- matrix(c

((t(df.Resid.M1.new$Q1)%*%df.Resid.M1.new$Q1)/72,

(t(df.Resid.M1.new$Q1)%*%df.Resid.M1.new$Q2)/72,

(t(df.Resid.M1.new$Q1)%*%df.Resid.M1.new$P1)/72,

(t(df.Resid.M1.new$Q1)%*%df.Resid.M1.new$P2)/72,

(t(df.Resid.M1.new$Q2)%*%df.Resid.M1.new$Q1)/72,

(t(df.Resid.M1.new$Q2)%*%df.Resid.M1.new$Q2)/72,

(t(df.Resid.M1.new$Q2)%*%df.Resid.M1.new$P1)/72,

(t(df.Resid.M1.new$Q2)%*%df.Resid.M1.new$P2)/72,

(t(df.Resid.M1.new$P1)%*%df.Resid.M1.new$Q1)/72,

(t(df.Resid.M1.new$P1)%*%df.Resid.M1.new$Q2)/72,

(t(df.Resid.M1.new$P1)%*%df.Resid.M1.new$P1)/72,

(t(df.Resid.M1.new$P1)%*%df.Resid.M1.new$P2)/72,

(t(df.Resid.M1.new$P2)%*%df.Resid.M1.new$Q1)/72,

(t(df.Resid.M1.new$P2)%*%df.Resid.M1.new$Q2)/72,

(t(df.Resid.M1.new$P2)%*%df.Resid.M1.new$P1)/72,

(t(df.Resid.M1.new$P2)%*%df.Resid.M1.new$P2)/72),

byrow=T, ncol=4)

colnames(new.vcov.M1) <- c("sigma_Q1Q1", "sigma_Q1Q2", "sigma_Q1P1", "sigma_Q1P2")

new.vcov.M1

new.vcov.M2 <- matrix(c

((t(df.Resid.M2.new$Q1)%*%df.Resid.M2.new$Q1)/72,

(t(df.Resid.M2.new$Q1)%*%df.Resid.M2.new$Q2)/72,

(t(df.Resid.M2.new$Q1)%*%df.Resid.M2.new$P1)/72,

(t(df.Resid.M2.new$Q1)%*%df.Resid.M2.new$P2)/72,

(t(df.Resid.M2.new$Q2)%*%df.Resid.M2.new$Q1)/72,

(t(df.Resid.M2.new$Q2)%*%df.Resid.M2.new$Q2)/72,

(t(df.Resid.M2.new$Q2)%*%df.Resid.M2.new$P1)/72,

(t(df.Resid.M2.new$Q2)%*%df.Resid.M2.new$P2)/72,

(t(df.Resid.M2.new$P1)%*%df.Resid.M2.new$Q1)/72,

(t(df.Resid.M2.new$P1)%*%df.Resid.M2.new$Q2)/72,

(t(df.Resid.M2.new$P1)%*%df.Resid.M2.new$P1)/72,

(t(df.Resid.M2.new$P1)%*%df.Resid.M2.new$P2)/72,

(t(df.Resid.M2.new$P2)%*%df.Resid.M2.new$Q1)/72,

(t(df.Resid.M2.new$P2)%*%df.Resid.M2.new$Q2)/72,

(t(df.Resid.M2.new$P2)%*%df.Resid.M2.new$P1)/72,

(t(df.Resid.M2.new$P2)%*%df.Resid.M2.new$P2)/72),

byrow=T, ncol=4)

colnames(new.vcov.M2) <- c("sigma_Q1Q1", "sigma_Q1Q2", "sigma_Q1P1", "sigma_Q1P2")

new.vcov.M2

library(MASS)

new.vcov.M1.I <- ginv(new.vcov.M1)

new.vcov.M2.I <- ginv(new.vcov.M2)

x.M1.M2=matrix(data=NA, nrow=72, ncol=1)

for(j in 1:72){

x.M1.M2[j,] <- ((Resid.Mat.M1[j,]%*%new.vcov.M1.I%*%Resid.Mat.M1.t[,j])

- (Resid.Mat.M2[j,]%*%new.vcov.M2.I%*%Resid.Mat.M2.t[,j]))^2

}

head(x.M1.M2)

normalizer.M1.M2 <- (sqrt(colSums(x.M1.M2)))/2

normalized.LR.M1.M2 <- 2*(-1735.88-(-1765.71))/normalizer.M1.M2

normalizer.M1.M2