library(spdep)
library(ggplot2)
library(spatialreg)

set.seed(123) 
S <- 500 # pocet simulacii
n <- 225 # pocet oblasti
lambda <- 0.6 # skutocna hodnota lambda

# skutocne hodnoty parametrov beta
beta0 <- 0
beta1 <- 1
beta2 <- 2
w <- cell2nb(sqrt(n), sqrt(n), type = 'queen') # vytvorenie matice vah
iw <- invIrM(w, lambda) # vypocita inverziu (I - rho*W)
beta0_hat <- vector(mode = "numeric", length = S) # vektor do ktoreho sa ulozia vysledky
beta1_hat <- vector(mode = "numeric", length = S) # vektor do ktoreho sa ulozia vysledky
beta2_hat <- vector(mode = "numeric", length = S) # vektor do ktoreho sa ulozia vysledky
lambda_hat <- vector(mode = "numeric", length = S) # vektor do ktoreho sa ulozia vysledky

# vytvorenie vysvetlujucich premennych x1 a x2
x1 <- 10 * rnorm(n, 1, 1) + seq(from = 1, to = 50, length = n)
x2 <- rnorm(n, 8, 1) + seq(from = 1, to = 7, length = n)
sigma = 1 * seq(from = 1, to = 15, length = n) # vektor simulujuci heteroskedasticitu
y <- x1 + 2 * x2 + iw %*% rnorm(n,mean=0,sd=sigma)

# grafy zavislosti y a x1, y a x2
df <- data.frame(y, x1, x2)
ggplot(df, aes(x = x1, y = y)) + geom_point() + theme_minimal()
ggplot(df, aes(x = x2, y = y)) + geom_point() + theme_minimal()


# v prvom riadku cyklu si definujeme skutocny priestorovy proces pomocou beta0, beta1, beta2, lambda
for (s in 1:S) {
  y <- x1 + 2 * x2 + iw %*% rnorm(n,mean=0,sd=sigma) # takto si definujeme skutocny priestorovy proces
  out <- errorsarlm(y ~ x1 + x2, listw = nb2listw(w), weights = 1 / sigma, method = 'LU') # odhad vazeneho SEM modelu
  beta0_hat[s] <- coef(out)[2] 
  beta1_hat[s] <- coef(out)[3] 
  beta2_hat[s] <- coef(out)[4] 
  lambda_hat[s] <- out$lambda 
}

residuals <- out$residuals

df <- data.frame(beta0_hat, beta1_hat, beta2_hat, lambda_hat)
df$beta0_hat <- beta0_hat
df$beta1_hat <- beta1_hat
df$beta2_hat <- beta2_hat

p <- ggplot(df, aes(x=beta0_hat)) + geom_density(color="darkblue", size = 1) + theme_minimal() + 
  xlab(expression(beta0_hat)) + geom_vline(xintercept = beta0, color = "red", size=1) 

q <- ggplot(df, aes(x=beta1_hat)) + geom_density(color="darkblue", size = 1) + theme_minimal() + 
  xlab(expression(beta1_hat)) + geom_vline(xintercept = beta1, color = "red", size=1) 

r <- ggplot(df, aes(x=beta2_hat)) + geom_density(color="darkblue", size = 1) + theme_minimal() + 
  xlab(expression(beta2_hat)) + geom_vline(xintercept = beta2, color = "red", size=1) 


s <- ggplot(df, aes(x=lambda_hat)) + geom_density(color="darkblue", size = 1) + theme_minimal() + 
  xlab(expression(hat(lambda))) + geom_vline(xintercept = lambda, color = "red", size=1) 

p
q
r
s

# tvar rezidui
df <- data.frame(residuals(out))
df$residuals <- residuals(out)
ggplot(df, aes(sample = residuals)) +
  stat_qq()  + theme_minimal() + 
  stat_qq_line()

