## ----setup, include=FALSE-----------------------------------------------------
rmarkdown::find_pandoc(version = '2.9.2.1')
knitr::opts_chunk$set(echo = TRUE, eval = TRUE)
knitr::knit_hooks$set(output = function(x, options) {
  x <- x
  paste0("\\begingroup\\small\n\\begin{verbatim}\n", x, "\\end{verbatim}\n\\endgroup\n")
})

## -----------------------------------------------------------------------------
library(QuantilePeer)
set.seed(123)  # Set seed for reproducibility
ngr  <- 50  # Number of subnets
nvec <- rep(30, ngr)  # Size of subnets
n    <- sum(nvec)

# Network matrix
G <- lapply(1:ngr, function(z) {
  Gz <- matrix(rbinom(nvec[z]^2, 1, 0.3), nvec[z], nvec[z])
  diag(Gz) <- 0
  # Adding isolated nodes (important for the structural model)
  niso <- sample(0:nvec[z], 1, prob = (nvec[z] + 1):1 / sum((nvec[z] + 1):1))
  if (niso > 0) {
    Gz[sample(1:nvec[z], niso), ] <- 0
  }
  Gz
})

X   <- cbind(rnorm(n), rpois(n, 2)); colnames(X) <- c("X1", "X2")

## -----------------------------------------------------------------------------
tau       <- seq(0, 1, 1/3) #quantile level
lambdatau <- c(0.1, 0.25, 0.2, 0.15) #lambda_tau
lambda2   <- 0.2 # lambda_2
beta      <- c(2, -0.5, 1)

# First dependent variable (reduced form without conformity)
y1        <- qpeer.sim(formula = ~ X, Glist = G, tau = tau, lambda = lambdatau,
                       beta = beta, structural = FALSE,
                       epsilon = rnorm(n, 0, 0.4))
y1        <- y1$y #qpeer.sim returns a list of several object including y

# Second dependent variable (structural form with conformity)
y2        <- qpeer.sim(formula = ~ X, Glist = G, tau = tau,
                       lambda = c(lambda2, lambdatau), beta = beta,
                       structural = TRUE, epsilon = rnorm(n, 0, 0.4))
y2        <- y2$y

## -----------------------------------------------------------------------------
# First instrument set
Z1  <- qpeer.inst(formula = ~ X, Glist = G, tau = seq(0, 1, 0.1),  
                  max.distance = 2, checkrank = TRUE) # finer subdivision
Z1  <- Z1$instruments #qpeer.inst returns a list of several object

# Second instrument set: y1 is used to order X
Z21 <- qpeer.inst(formula = y1 ~ X, Glist = G, tau = tau, max.distance = 2,
                  checkrank = TRUE)
qy1 <- Z21$qy #quantile of y among peers                 
Z21 <- Z21$instruments

# Second instrument set: y2 is used to order X
Z22 <- qpeer.inst(formula = y2 ~ X, Glist = G, tau = tau, max.distance = 2,
                  checkrank = TRUE)
qy2 <- Z22$qy #quantile of y among peers                 
Z22 <- Z22$instruments

## -----------------------------------------------------------------------------
QtR1 <- qpeer(formula = y1 ~ X, excluded.instruments = ~ Z1, Glist = G,
              tau = tau)
summary(QtR1, diagnostic = TRUE)

QtR2 <- qpeer(formula = y1 ~ X, excluded.instruments = ~ Z1 + Z21, Glist = G,
              tau = tau)
summary(QtR2)

QtR3 <- qpeer(formula = y2 ~ X, excluded.instruments = ~ Z1, Glist = G,
              tau = tau)
summary(QtR3)

QtR4 <- qpeer(formula = y2 ~ X, excluded.instruments = ~ Z1 + Z21, Glist = G,
             tau = tau)
summary(QtR4)

## -----------------------------------------------------------------------------
QtS1 <- qpeer(formula = y1 ~ X, excluded.instruments = ~ Z1, Glist = G,
              tau = tau, structural = TRUE)
summary(QtS1, diagnostic = TRUE)

QtS2 <- qpeer(formula = y1 ~ X, excluded.instruments = ~ Z1 + Z21, Glist = G,
              tau = tau, structural = TRUE)
summary(QtS2)

QtS3 <- qpeer(formula = y2 ~ X, excluded.instruments = ~ Z1, Glist = G,
              tau = tau, structural = TRUE)
summary(QtS3)

QtS4 <- qpeer(formula = y2 ~ X, excluded.instruments = ~ Z1 + Z22, Glist = G,
              tau = tau, structural = TRUE)
summary(QtS4)

## -----------------------------------------------------------------------------
QtR5  <- qpeer(formula = y1 ~ X, excluded.instruments = ~ Z1, Glist = G,
               tau = tau, structural = FALSE, estimator = "gmm.optimal",
               HAC = "cluster", fixed.effects = "separate")
summary(QtR5)

QtS5 <- qpeer(formula = y2 ~ X, excluded.instruments = ~ Z1, Glist = G,
              tau = tau, structural = TRUE, estimator = "gmm.optimal",
              HAC = "cluster", fixed.effects = "separate")
summary(QtS5)

## -----------------------------------------------------------------------------
qpeer.test(QtR1, which = "uniform")
qpeer.test(QtS5, which = "decreasing")

## -----------------------------------------------------------------------------
qpeer.test(QtR1, QtR2, which = "wald")
qpeer.test(QtS3, QtS4, which = "sargan")

## -----------------------------------------------------------------------------
# Estimating QtS6 with a misspecified tau
QtS6 <- qpeer(formula = y2 ~ X, excluded.instruments = ~ Z1, Glist = G,
              tau = c(0, 1), structural = TRUE, estimator = "gmm.optimal",
              HAC = "cluster", fixed.effects = "separate")
qpeer.test(model1 = QtS6, model2 = QtS5, which = "encompassing")

## -----------------------------------------------------------------------------
qpeer.test(model1 = QtS5, model2 = QtS6, which = "encompassing")

## -----------------------------------------------------------------------------
ngr  <- 50  # Number of subnets
nvec <- rep(30, ngr)  # Size of subnets
n    <- sum(nvec)

# Network matrix
G <- lapply(1:ngr, function(z) {
  Gz <- matrix(rbinom(nvec[z]^2, 1, 0.3), nvec[z], nvec[z])
  diag(Gz) <- 0
  # Adding isolated nodes (important for the structural model)
  niso <- sample(0:nvec[z], 1, prob = (nvec[z] + 1):1 / sum((nvec[z] + 1):1))
  if (niso > 0) {
    Gz[sample(1:nvec[z], niso), ] <- 0
  }
  Gz
  rs   <- rowSums(Gz); rs[rs == 0] <- 1
  Gz   <- Gz/rs # rowSums are normalized to one
})

X   <- cbind(rnorm(n), rpois(n, 2)); colnames(X) <- c("X1", "X2")

lambda  <- 0.55
lambda2 <- 0.2
beta    <- c(2.5, -0.5, 1)
rho     <- -3

# First dependent variable (reduced form without conformity)
y1      <- cespeer.sim(formula = ~ X, Glist = G, rho = rho,
                       lambda = lambda, beta = beta, structural = FALSE,
                       epsilon = rnorm(n, 0, 0.4))
y1      <- y1$y #cespeer.sim returns a list of several object including y

# Second dependent variable (structural form with conformity)
y2      <-  cespeer.sim(formula = ~ X, Glist = G, rho = rho,
                        lambda = c(lambda2, lambda), beta = beta,
                        structural = TRUE, epsilon = rnorm(n, 0, 0.4))
y2      <- y2$y

## -----------------------------------------------------------------------------
yhat1 <- fitted(lm(y1 ~ X))
yhat2 <- fitted(lm(y2 ~ X))

## -----------------------------------------------------------------------------
cesR1 <- cespeer(formula = y1 ~ X, instrument = ~ yhat1, Glist = G,
                 structural = FALSE)
summary(cesR1)
cesS1 <- cespeer(formula = y2 ~ X, instrument = ~ yhat2, Glist = G,
                 structural = TRUE)
summary(cesS1)

## -----------------------------------------------------------------------------
Z    <- qpeer.inst(formula = ~ X, Glist = G, tau = seq(0, 1, 0.1),  
                  max.distance = 2, checkrank = TRUE)$instruments
QtR1 <- qpeer(formula = y1 ~ X, excluded.instruments = ~ Z, Glist = G,
              tau = seq(0, 1, 1/3))
qpeer.test(QtR1, which = "decreasing")

## -----------------------------------------------------------------------------
library(PartialNetwork)
GX  <- peer.avg(G, X)
GGX <- peer.avg(G, GX)
Ot1 <- linpeer(formula = y2 ~ X, excluded.instruments = ~ GX + GGX, Glist = G,
               structural = TRUE, estimator = "gmm.optimal",
               fixed.effects = "separate", HAC = "cluster")
summary(Ot1, diagnostic = TRUE)

## -----------------------------------------------------------------------------
Gy2 <- peer.avg(G, y2)
qy2 <- qpeer.inst(formula = y2 ~ 1, Glist = G, tau = c(0, 1))$qy #min and max
Ot2 <- genpeer(formula = y2 ~ X, excluded.instruments = ~ Z + GX + GGX,
               endogenous.variables = ~ Gy2 + qy2, Glist = G, structural = TRUE,
               estimator = "gmm.optimal", fixed.effects = "separate",
               HAC = "cluster") # includes average, min, and max of peers
summary(Ot2, diagnostic = TRUE)