plot.random.effects <- function(fit, h, mode, distribution, name){
  
  p <- length(fit$theta.hat)
  N.stations <- N.stations + 1*(mode == "dyadic" & distribution != "poisson")
  u.hat <- tibble(u.out.hat = fit$theta.hat[(p - 2*(N.stations) + 1):(p - N.stations)],
                  u.in.hat = fit$theta.hat[(p - N.stations + 1):(p)])
  
  g <- ggplot(u.hat, aes(x = u.out.hat, y = u.in.hat)) + geom_point() +
    labs(title = paste0("Fitted Random Effects for h = ", h), x = expression(widehat(u)^out),
         y = expression(widehat(u)^"in"))
  
  pdf(file = paste0("Plots/Model/random_effects_", name, ".pdf"), width = 6, height = 4)
  print(g)
  dev.off()
  setEPS()
  postscript(paste0("Plots/Model/random_effects_", name, ".eps"), width = 6, height = 4)
  print(g)
  dev.off()
  
}

plot.smooth.effects <- function(fit, design, name){
  
  for (m in 1:length(smooths)) {
    z <- switch(smooths[m],
               "seas" = design$z[[m]],
               seq(min(design$z[[m]]), max(design$z[[m]]), length.out = 100))
    B <- switch(smooths[m], 
               "seas" = design$B[[m]],
               splineDesign(knots = design$k[[m]], x = z, ord = design$l[m] + 1, outer.ok = TRUE))
    if (smooths[m] != "seas"){
      B <- sweep(B, 2, colMeans(B))[, -1]
    }
    gamma.hat <- fit$theta.hat[design$ind.gamma[[m]]]
    s.hat <- as.vector(B%*%gamma.hat)
    conf.band <- confidence.band(gamma = gamma.hat, cov = as.matrix(fit$V.hat[design$ind.gamma[[m]], design$ind.gamma[[m]]]), 
                                z = z, B = B)
    
    values <- tibble(z = z, val = s.hat, lower = conf.band$lower, upper = conf.band$upper)
    
    g <- ggplot(values, aes(x = z, y = val)) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.3) + geom_line(color = "red") + 
      labs(x = bquote(.(smooths[m])[t]), y = bquote(s[.(m)] ~ (.(smooths[m])[t])), 
           title = paste0("Estimated Smooth Effect of ", smooths[m]))
    if (smooths[m] == "seas"){
      g <- g + scale_x_continuous(breaks = seq(1/24, 23/24, 1/12), 
                                  label = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
    }

    pdf(file = paste0("Plots/Model/smooths_", name, "_", smooths[m], ".pdf"), width = 6, height = 4)
    print(g)
    dev.off()
  }
}

plot.cum.by.station <- function(observed, estimated, design, lim, name){
  
  counts <- tibble(obs = c(colSums(observed$outgoing[design$H, ]), colSums(observed$incoming[design$H, ])),
                      est = c(colSums(estimated$N.out.hat[, 1:N.stations]), colSums(estimated$N.in.hat[, 1:N.stations])),
                      kind = c(rep("outgoing", N.stations), rep("incoming", N.stations)))
  counts$kind <- factor(counts$kind, levels = c("outgoing", "incoming"))
  
  g <- ggplot(counts, aes(x = obs, y = est)) + geom_point(aes(color = kind, shape = kind)) + 
    geom_abline(intercept = 0, slope = 1) +
    scale_y_continuous(breaks = seq(0, lim, 500), limits = c(0, lim)) +
    scale_x_continuous(breaks = seq(0, lim, 500), limits = c(0, lim)) +
    theme(legend.justification = c(0.01, 0.99), legend.position = c(0.01, 0.99), legend.title = element_blank()) + 
    labs(x = "Cumulated Observed Trip Count", y = "Cumulated Estimated Trip Count", title = "Cumulated Trip Counts by Station") 
  
  pdf(file = paste0("Plots/Model/by_station_", name, ".pdf"), width = 6, height = 4)
  print(g)
  dev.off()
  setEPS()
  postscript(paste0("Plots/Model/by_station_", name, ".eps"), width = 6, height = 4)
  print(g)
  dev.off()
}

plot.cum.by.day <- function(observed, estimated, design, lim, name){
  
  counts <- tibble(index = rep(1:length(design$H), 2), 
                      counts = c(rowSums(observed$outgoing[design$H, ]), 
                                 rowSums(estimated$N.out.hat[, 1:N.stations])),
                      kind = c(rep("observed", length(design$H)), 
                               rep("estimated", length(design$H))))
  
  g <- ggplot(counts, aes(x = index, y = counts)) + geom_line(aes(linetype = kind, color = kind)) + 
    theme(legend.justification = c(0.01, 0.99), legend.position = c(0.01, 0.99), legend.title = element_blank()) + 
    labs(x = "Day Index", y = "Cumulated Trip Count", title = "Cumulated Trip Counts by Day Index") + 
    scale_y_continuous(breaks = seq(0, lim, 100), limits = c(0, lim)) +
    scale_x_continuous(breaks = seq(0, length(design$H), 25))
  
  pdf(file = paste0("Plots/Model/by_day_", name, ".pdf"), width = 6, height = 4)
  print(g)
  dev.off()
  setEPS()
  postscript(paste0("Plots/Model/by_day_", name, ".eps"), width = 6, height = 4)
  print(g)
  dev.off()
}

plot.residuals <- function(observed, estimated, design, name){
  
  pearson.out <- pearson.in <- matrix(NA, length(design$H), N.stations)
  for (i in 1:N.stations) {
    ind <- which(design$H >= stations$ActiveFrom[i]) 
    pearson.out[ind, i] <- (observed$outgoing[design$H[ind], i] - estimated$N.out.hat[ind, i])/
      (sqrt(estimated$N.out.hat[ind, i]))
    pearson.in[ind, i] <- (observed$outgoing[design$H[ind], i] - estimated$N.in.hat[ind, i])/
      (sqrt(estimated$N.in.hat[ind, i]))
  }  
  pearson.out <- pearson.out[which(!is.na(pearson.out))]
  pearson.in <- pearson.in[which(!is.na(pearson.in))]
  
  residuals <- data.frame(resid = c(pearson.out, pearson.in), 
                      kind = c(rep("outgoing", length(pearson.out)), 
                               rep("incoming", length(pearson.in))))
  residuals$kind <- factor(residuals$kind, levels = c("outgoing", "incoming"))
  
  g <- ggplot(residuals, aes(x = resid)) + 
    geom_density(aes(color = kind, linetype = kind, fill = kind), alpha = 0.1) + 
    geom_hline(color = "white", yintercept = 0) + 
    scale_x_continuous(limits = c(-3, 4), breaks = seq(-3, 4, 1)) + 
    scale_y_continuous(breaks = seq(0, 1, 0.1)) + 
    theme(legend.justification = c(0.99, 0.99), legend.position = c(0.99, 0.99), legend.title = element_blank()) + 
    labs(x = "Pearson Residual", y = "Density", "Density of Pearson Residuals")
  
  pdf(file = paste0("Plots/Model/residuals_", name, ".pdf"), width = 6, height = 4)
  print(g)
  dev.off()
}

plot.distributions <- function() {
  
  files <- list.files(path = "Fits/")
  files <- files[intersect(grep("_corrected", files), grep("_station", files))]
  hour <- as.numeric(sub("\\_.*", "", files))
  skellam <- grep("_skellam", files)
  tskellam <- grep("_tskellam", files)
  pois <- grep("_poisson", files)
  distributions <- rep("", length(files))
  distributions[skellam] <- "skellam"
  distributions[tskellam] <- "tskellam"
  distributions[pois] <- "poisson"
  
  Delta.day.out <- Delta.station.out <- vector("list", length(files))
  for (i in 1:length(files)) {
    design <- readRDS(paste0("Designs/", files[i]))
    fit <- readRDS(paste0("Fits/", files[i]))
    observed <-  switch(interventions[i],
                        "none" = N$outgoing[design$H, ],
                        N.wo.P$outgoing[design$H, ])
    estimated <- estimate.flows(fit = fit, design = design, mode = "station")$N.out.hat
    
    Delta.day.out[[i]] <- sum(abs((rowSums(observed - estimated)))/sum(observed))
    Delta.station.out[[i]] <- sum(abs((colSums(observed - estimated)))/sum(observed))
  }
  
  range.delta <- range(c(unlist(Delta.day.out), unlist(Delta.station.out)))
  df.day <- tibble(h = factor(hour), 
                   delta = unlist(Delta.day.out, use.names = FALSE),
                   distribution = factor(distributions, levels = c("skellam", "tskellam", "poisson")))
  df.station <- tibble(h = factor(hour), 
                       delta = unlist(Delta.station.out, use.names = FALSE),
                       distribution = factor(distributions, levels = c("skellam", "tskellam", "poisson")))
  
  g <- ggplot() +
    geom_point(data = df.day, aes(x = h, y = delta, color = distribution, shape = distribution), size = 2.3) +
    geom_line(data = df.day, aes(x = h, y = delta, group = distribution, color = distribution)) + 
    labs(x = "h", y = expression(Delta[time]^out*(h)), title = "Time Specific Estimation Error", 
         color = "Model", shape = "Model") +
    scale_y_continuous(limits = c(0, max(range.delta))) +
    theme(legend.justification = c(0.01, 0.99), legend.position = c(0.01, 0.99),
          legend.direction = "horizontal") + 
    scale_color_hue(labels = c("Skellam", "T-Skellam", "Poisson")) + 
    scale_shape(labels = c("Skellam", "T-Skellam", "Poisson"))
  
  pdf(file = "Plots/Model/delta_day_out_distributions.pdf", width = 6, height = 4)
  print(g)
  dev.off()
  setEPS()
  postscript("Plots/Model/delta_day_out_distributions.eps", width = 6, height = 4)
  print(g)
  dev.off()
  
  g <- ggplot() +
    geom_point(data = df.station, aes(x = h, y = delta, color = distribution, shape = distribution), size = 2.3) +
    geom_line(data = df.station, aes(x = h, y = delta, group = distribution, color = distribution)) + 
    labs(x = "h", y = expression(Delta[stat]^out*(h)), title = "Station Specific Estimation Error", 
         color = "Model", shape = "Model") +
    scale_y_continuous(limits = c(0, max(range.delta))) +
    theme(legend.justification = c(0.01, 0.99), legend.position = c(0.01, 0.99),
          legend.direction = "horizontal") + 
    scale_color_hue(labels = c("Skellam", "T-Skellam", "Poisson")) + 
    scale_shape(labels = c("Skellam", "T-Skellam", "Poisson"))
  
  pdf(file = "Plots/Model/delta_station_out_distributions.pdf", width = 6, height = 4)
  print(g)
  dev.off()
  setEPS()
  postscript("Plots/Model/delta_station_out_distributions.eps", width = 6, height = 4)
  print(g)
  dev.off()
}

plot.interventions <- function() {
  
  files <- list.files(path = "Fits/")
  files <- files[intersect(grep("_skellam", files), grep("_station", files))]
  hour <- as.numeric(sub("\\_.*", "", files))
  corrected <- grep("corrected_", files)
  threshold <- grep("threshold", files)
  none <- grep("none", files)
  interventions <- rep("", length(files))
  interventions[corrected] <- "corrected"
  interventions[threshold] <- "threshold"
  interventions[none] <- "none"
  
  Delta.day.out <- Delta.station.out <- vector("list", length(files))
  for (i in 1:length(files)) {
    design <- readRDS(paste0("Designs/", files[i]))
    fit <- readRDS(paste0("Fits/", files[i]))
    observed <-  switch(interventions[i],
                        "none" = N$outgoing[design$H, ],
                        N.wo.P$outgoing[design$H, ])
    estimated <- estimate.flows(fit = fit, design = design, mode = "station")$N.out.hat
    
    Delta.day.out[[i]] <- sum(abs((rowSums(observed - estimated)))/sum(observed))
    Delta.station.out[[i]] <- sum(abs((colSums(observed - estimated)))/sum(observed))
  }
  
  range.delta <- range(c(unlist(Delta.day.out), unlist(Delta.station.out)))
  df.day <- tibble(h = factor(hour), 
                   delta = unlist(Delta.day.out, use.names = FALSE),
                   interventions = factor(interventions, levels = c("corrected", "threshold", "none")))
  df.station <- tibble(h = factor(hour), 
                       delta = unlist(Delta.station.out, use.names = FALSE),
                       interventions = factor(interventions, levels = c("corrected", "threshold", "none")))
  
  P.h <- rep(0, 12)
  for (h in 9:20) {
    H <- readRDS(paste0("Designs/", files[min(which(hour == h))]))$H
    P.h[h - 8] <- sum(P$outgoing[H, ])/sum(N$outgoing[H, ])
  }
  df.interventions <- tibble(h = factor(9:20), share.P = P.h)
  
  g <- ggplot() + geom_col(data = df.interventions, aes(x = h, y = share.P*4, fill = "Share of \nProvider Interventions"), alpha = 0.3) + 
    scale_fill_manual(values = "black") +
    geom_point(data = df.day, aes(x = h, y = delta, color = interventions, shape = interventions), size = 2.3) +
    geom_line(data = df.day, aes(x = h, y = delta, group = interventions, color = interventions)) + 
    labs(x = "h", y = expression(Delta[time]^out*(h)), title = "Time Specific Estimation Error", 
         color = "Treatment of \nProvider Interventions", shape = "Treatment of \nProvider Interventions", fill = "") +
    scale_y_continuous(limits = c(0, max(range.delta)), 
                       sec.axis = sec_axis(~ ./4, name = "Share of Provider Interventions", breaks = seq(0, 0.25, 0.05))) +
    theme(legend.justification = c(0.99, 0.99), legend.position = c(0.99, 0.99)) + 
    guides(color = guide_legend(order = 1), shape = guide_legend(order = 1), fill = guide_legend(order = 2))
  
  pdf(file = "Plots/Model/delta_day_out_interventions.pdf", width = 6, height = 4)
  print(g)
  dev.off()
  
  g <- ggplot() + geom_col(data = df.interventions, aes(x = h, y = share.P*4, fill = "Share of \nProvider Interventions"), alpha = 0.3) + 
    scale_fill_manual(values = "black") +
    geom_point(data = df.station, aes(x = h, y = delta, color = interventions, shape = interventions), size = 2.3) +
    geom_line(data = df.station, aes(x = h, y = delta, group = interventions, color = interventions)) + 
    labs(x = "h", y = expression(Delta[stat]^out*(h)), title = "Station Specific Estimation Error", 
         color = "Treatment of \nProvider Interventions", shape = "Treatment of \nProvider Interventions", fill = "") +
    scale_y_continuous(limits = c(0, max(range.delta)), 
                       sec.axis = sec_axis(~ ./4, name = "Share of Provider Interventions", breaks = seq(0, 0.25, 0.05))) +
    theme(legend.justification = c(0.99, 0.99), legend.position = c(0.99, 0.99)) + 
    guides(color = guide_legend(order = 1), shape = guide_legend(order = 1), fill = guide_legend(order = 2))
  
  pdf(file = "Plots/Model/delta_station_out_interventions.pdf", width = 6, height = 4)
  print(g)
  dev.off()
}

plot.cum.dist <- function(observed, estimations.skellam, estimations.poisson, design){
  # counts by connection
  by.connection.observed <- apply(observed$dyadic[design$H, , ], 2:3, sum)
  by.connection.skellam <- apply(estimations.skellam, 2:3, sum) 
  by.connection.poisson <- apply(estimations.poisson, 2:3, sum)

  # counts by distance between stations
  bins.distance <- seq(0, max(distances), 0.1)
  N.distance.observed <- rep(0, length(bins.distance))
  N.distance.poisson <- N.distance.skellam <- rep(0, length(bins.distance))
  for (k in 1:N.stations) {
    for (l in 1:N.stations) {
      dist <- ceiling(distances[k, l]*10)
      N.distance.observed[dist + 1] = N.distance.observed[dist + 1] + by.connection.observed[k, l]
      N.distance.skellam[dist + 1] = N.distance.skellam[dist+ 1] + by.connection.skellam[k, l]
      N.distance.poisson[dist + 1] = N.distance.poisson[dist + 1] + by.connection.poisson[k, l]
    }
  }
  
  # cumulative distribution function
  cum.dist.observed <- cumsum(N.distance.observed[2:length(N.distance.observed)])/sum(N.distance.observed[2:length(N.distance.observed)])
  cum.dist.skellam <- cumsum(N.distance.skellam[2:length(N.distance.skellam)])/sum(N.distance.skellam[2:length(N.distance.skellam)])
  cum.dist.poisson <- cumsum(N.distance.poisson[2:length(N.distance.poisson)])/sum(N.distance.poisson[2:length(N.distance.poisson)])
  
  r <- length(bins.distance[2:length(bins.distance)])
  df <- tibble(x = rep(bins.distance[2:length(bins.distance)], 3),
                   y = c(cum.dist.observed, cum.dist.skellam, cum.dist.poisson),
                   kind = c(rep("Observed", r), rep("Skellam", r), rep("Poisson", r)))
  df$kind <- factor(df$kind, levels = c("Observed", "Skellam", "Poisson"))
  
  g <- ggplot(df, aes(x = x, y = y)) + geom_line(aes(color = kind, linetype = kind)) + 
    labs(x = "d", y = expression(P(D <= d)), color = "", linetype = "") + 
    scale_x_continuous(breaks = seq(0, 12, 2)) + 
    theme(legend.justification = c(0.01, 0.99), legend.position = c(0.01, 0.99))
  
  pdf(file = "Plots/Model/cum_dist.pdf", width = 6, height = 4)
  print(g)
  dev.off()
  setEPS()
  postscript("Plots/Model/cum_dist.eps", width = 6, height = 4)
  print(g)
  dev.off()
}

plot.pois.vs.skellam <- function(estimations.poisson, estimations.skellam, design){
  # sum over all days
  by.connection.poisson <- apply(estimations.poisson, 2:3, sum)
  by.connection.skellam <- apply(estimations.skellam, 2:3, sum) 
  
  df <- tibble(x = as.vector(by.connection.poisson), y = as.vector(by.connection.skellam))
  lim <- max(by.connection.poisson, by.connection.skellam)
  
  g <- ggplot(df, aes(x = x, y = y)) + geom_point() + 
    scale_x_continuous(limits = c(0, lim)) +
    scale_y_continuous(limits = c(0, lim)) +
    labs(x = expression(hat(mu)["ij,"*bold(.)]^Poisson), y = expression(hat(mu)["ij,"*bold(.)]^Skellam)) + 
    geom_abline(slope = 1, intercept = 0)
  
  pdf(file = "Plots/Model/pois_vs_skellam.pdf", width = 6, height = 4)
  print(g)
  dev.off()
  setEPS()
  postscript("Plots/Model/pois_vs_skellam.eps", width = 6, height = 4)
  print(g)
  dev.off()
}

plot.check.poisson <- function(observed, estimated, design, name){
  prob.0 <- prob.1 <- prob.2 <- rep(0, length(design$H))
  prob.0.hat <- prob.1.hat <- prob.2.hat <- rep(0, length(design$H))
  for (d in 1:length(design$H)) {
    prob.0[d] <- length(which(as.vector(N$dyadic[design$H[d], , ]) == 0))/(N.stations^2)
    prob.1[d] <- length(which(as.vector(N$dyadic[design$H[d], , ]) == 1))/(N.stations^2)
    prob.2[d] <- length(which(as.vector(N$dyadic[design$H[d], , ]) >= 2))/(N.stations^2)
    prob.0.hat[d] <- mean(exp(-estimated[d, , ]))
    prob.1.hat[d] <- mean(estimated[d, , ]*exp(-estimated[d, , ]))
    prob.2.hat[d] <- mean(1 - exp(-estimated[d, , ]) - estimated[d, , ]*exp(-estimated[d, , ]))
  }
  
  df.0 <- tibble(x = prob.0, y = prob.0.hat)
  df.1 <- tibble(x = prob.1, y = prob.1.hat)
  df.2 <- tibble(x = prob.2, y = prob.2.hat)
  
  g.0 <- ggplot(df.0, aes(x = x, y = y)) + geom_point() + geom_abline(slope = 1, intercept = 0) + 
    labs(y = expression(widehat(P)(N[t] == 0)), x = expression("Share of "*N[t] == 0)) + 
    scale_x_continuous(limits = c(min(prob.0, prob.0.hat), 1)) + 
    scale_y_continuous(limits = c(min(prob.0, prob.0.hat), 1)) + 
    theme(axis.text = element_text(size = 12), axis.title = element_text(size = 12))
  g.1 <- ggplot(df.1, aes(x = x, y = y)) + geom_point() + geom_abline(slope = 1, intercept = 0) + 
    labs(y = expression(widehat(P)(N[t] == 1)), x = expression("Share of "*N[t] == 1)) + 
    scale_x_continuous(limits = c(0, max(prob.1, prob.1.hat))) + 
    scale_y_continuous(limits = c(0, max(prob.1, prob.1.hat))) + 
    theme(axis.text = element_text(size = 12), axis.title = element_text(size = 12))
  g.2 <- ggplot(df.2, aes(x = x, y = y)) + geom_point() + geom_abline(slope = 1, intercept = 0) + 
    labs(y = expression(widehat(P)(N[t] >= 2)), x = expression("Share of "*N[t] >= 2)) +  
    scale_x_continuous(limits = c(0, max(prob.2, prob.2.hat))) +  
    scale_y_continuous(limits = c(0, max(prob.2, prob.2.hat))) + 
    theme(axis.text = element_text(size = 12), axis.title = element_text(size = 12))
  
  pdf(file = paste0("Plots/Model/probabilities", name, ".pdf"), width = 12, height = 8)
  print(grid.arrange(g.0, g.1, g.2, ncol = 3))
  dev.off()
  setEPS()
  postscript(paste0("Plots/Model/probabilities", name, ".eps"), width = 12, height = 8)
  print(grid.arrange(g.0, g.1, g.2, ncol = 3))
  dev.off()
}

plot.simulation <- function(results, S, N.stations) {
  
  beta.0.hat <- beta.1.hat <- beta.2.hat <- Sigma.11.hat <- Sigma.22.hat <- Sigma.12.hat <-matrix(0, S, length(results))
  for (b in 1:length(results)) {
    for (s in 1:S) {
      beta.0.hat[s, b] = results[[b]][[s]]$theta.hat[1]
      beta.1.hat[s, b] = results[[b]][[s]]$theta.hat[2]
      beta.2.hat[s, b] = results[[b]][[s]]$theta.hat[3]
      Sigma.11.hat[s, b] = results[[b]][[s]]$Sigma.star[1, 1]
      Sigma.22.hat[s, b] = results[[b]][[s]]$Sigma.star[2, 2]
      Sigma.12.hat[s, b] = results[[b]][[s]]$Sigma.star[1, 2]
    }
  }
  
  pdf(file = paste0("Plots/Simulation/BoxplotParameters_N", N.stations, "_S", S, ".pdf"), width = 11/1.5, height = 6/1.5)
  par(mar=c(2.5, 4.5, 0.5, 1))
  par(mfrow = c(2, 3))
  boxplot(beta.0.hat, names = c(expression(widehat(beta)[0]^(1)), 
                                expression(widehat(beta)[0]^(2)), 
                                expression(widehat(beta)[0]^(3)), 
                                expression(widehat(beta)[0]^(4))))
  lines(c(0.5, 1.5), c(-5, -5), lty = 2, col = "red", lwd = 2)
  lines(c(1.5, 2.5), c(-5, -5), lty = 2, col = "red", lwd = 2)
  lines(c(2.5, 3.5), c(-5, -5), lty = 2, col = "red", lwd = 2)
  lines(c(3.5, 4.5), c(-5, -5), lty = 2, col = "red", lwd = 2)
  boxplot(beta.1.hat, names = c(expression(widehat(beta)[1]^(1)), 
                                expression(widehat(beta)[1]^(2)), 
                                expression(widehat(beta)[1]^(3)), 
                                expression(widehat(beta)[1]^(4))))
  lines(c(0.5, 1.5), c(1, 1), lty = 2, col = "red", lwd = 2)
  lines(c(1.5, 2.5), c(-1, -1), lty = 2, col = "red", lwd = 2)
  lines(c(2.5, 3.5), c(1, 1), lty = 2, col = "red", lwd = 2)
  lines(c(3.5, 4.5), c(-1, -1), lty = 2, col = "red", lwd = 2)
  boxplot(beta.2.hat, names = c(expression(widehat(beta)[2]^(1)), 
                                expression(widehat(beta)[2]^(2)), 
                                expression(widehat(beta)[2]^(3)), 
                                expression(widehat(beta)[2]^(4))))
  lines(c(0.5, 1.5), c(-1, -1), lty = 2, col = "red", lwd = 2)
  lines(c(1.5, 2.5), c(1, 1), lty = 2, col = "red", lwd = 2)
  lines(c(2.5, 3.5), c(1, 1), lty = 2, col = "red", lwd = 2)
  lines(c(3.5, 4.5), c(-1, -1), lty = 2, col = "red", lwd = 2)
  boxplot(Sigma.11.hat, names = c(expression(widehat(sigma)[1]^2), 
                                  expression(widehat(sigma)[1]^2), 
                                  expression(widehat(sigma)[1]^2), 
                                  expression(widehat(sigma)[1]^2)))
  lines(c(0, 5), c(1, 1), lty = 2, col = "red", lwd = 2)
  boxplot(Sigma.22.hat, names = c(expression(widehat(sigma)[2]^2), 
                                  expression(widehat(sigma)[2]^2), 
                                  expression(widehat(sigma)[2]^2), 
                                  expression(widehat(sigma)[2]^2)))
  lines(c(0, 5), c(1, 1), lty = 2, col = "red", lwd = 2)
  boxplot(Sigma.12.hat, names = c(expression(widehat(sigma)[12]^2), 
                                  expression(widehat(sigma)[12]^2), 
                                  expression(widehat(sigma)[12]^2), 
                                  expression(widehat(sigma)[12]^2)))
  lines(c(0, 5), c(0.9, 0.9), lty = 2, col = "red", lwd = 2)
  dev.off()
  
  setEPS()
  postscript(paste0("Plots/Simulation/BoxplotParameters_N", N.stations, "_S", S, ".eps"), width = 11/1.5, height = 6/1.5)
  par(mar=c(2.5, 4.5, 0.5, 1))
  par(mfrow = c(2, 3))
  boxplot(beta.0.hat, names = c(expression(widehat(beta)[0]^(1)), 
                                expression(widehat(beta)[0]^(2)), 
                                expression(widehat(beta)[0]^(3)), 
                                expression(widehat(beta)[0]^(4))))
  lines(c(0.5, 1.5), c(-5, -5), lty = 2, col = "red", lwd = 2)
  lines(c(1.5, 2.5), c(-5, -5), lty = 2, col = "red", lwd = 2)
  lines(c(2.5, 3.5), c(-5, -5), lty = 2, col = "red", lwd = 2)
  lines(c(3.5, 4.5), c(-5, -5), lty = 2, col = "red", lwd = 2)
  boxplot(beta.1.hat, names = c(expression(widehat(beta)[1]^(1)), 
                                expression(widehat(beta)[1]^(2)), 
                                expression(widehat(beta)[1]^(3)), 
                                expression(widehat(beta)[1]^(4))))
  lines(c(0.5, 1.5), c(1, 1), lty = 2, col = "red", lwd = 2)
  lines(c(1.5, 2.5), c(-1, -1), lty = 2, col = "red", lwd = 2)
  lines(c(2.5, 3.5), c(1, 1), lty = 2, col = "red", lwd = 2)
  lines(c(3.5, 4.5), c(-1, -1), lty = 2, col = "red", lwd = 2)
  boxplot(beta.2.hat, names = c(expression(widehat(beta)[2]^(1)), 
                                expression(widehat(beta)[2]^(2)), 
                                expression(widehat(beta)[2]^(3)), 
                                expression(widehat(beta)[2]^(4))))
  lines(c(0.5, 1.5), c(-1, -1), lty = 2, col = "red", lwd = 2)
  lines(c(1.5, 2.5), c(1, 1), lty = 2, col = "red", lwd = 2)
  lines(c(2.5, 3.5), c(1, 1), lty = 2, col = "red", lwd = 2)
  lines(c(3.5, 4.5), c(-1, -1), lty = 2, col = "red", lwd = 2)
  boxplot(Sigma.11.hat, names = c(expression(widehat(sigma)[1]^2), 
                                  expression(widehat(sigma)[1]^2), 
                                  expression(widehat(sigma)[1]^2), 
                                  expression(widehat(sigma)[1]^2)))
  lines(c(0, 5), c(1, 1), lty = 2, col = "red", lwd = 2)
  boxplot(Sigma.22.hat, names = c(expression(widehat(sigma)[2]^2), 
                                  expression(widehat(sigma)[2]^2), 
                                  expression(widehat(sigma)[2]^2), 
                                  expression(widehat(sigma)[2]^2)))
  lines(c(0, 5), c(1, 1), lty = 2, col = "red", lwd = 2)
  boxplot(Sigma.12.hat, names = c(expression(widehat(sigma)[12]^2), 
                                  expression(widehat(sigma)[12]^2), 
                                  expression(widehat(sigma)[12]^2), 
                                  expression(widehat(sigma)[12]^2)))
  lines(c(0, 5), c(0.9, 0.9), lty = 2, col = "red", lwd = 2)
  dev.off()
}