# Aula 1 ------------------------------------------------------------------
# Exemplo de cálculo de matriz de Covariâncias e de Correlações amostrais
# Exemplo 1 - cov + corr + MahalanobisArquivo
rm(list = ls()) #limpando a memória do R

# Definindo os vetores 
y1 <- c(72,60,56,41,32,30,39,42,37,33,32,63,54,47,91,56,79,81,78,46,39,32,60,35,39,50,43,48)
y2 <- c(66,53,57,29,32,35,39,43,40,29,30,45,46,51,79,68,65,80,55,38,35,30,50,37,36,34,37,54)
y3 <- c(76,66,64,36,35,34,31,31,31,27,34,74,60,52,100,47,70,68,67,37,34,30,67,48,39,37,39,57)
y4 <- c(77,63,58,38,36,26,27,25,25,36,28,63,52,43,75,50,61,58,60,38,37,32,54,39,31,40,50,43)

#Organizando os dados
Y <- cbind(y1,y2,y3,y4) #Criando uma matriz com os vetores de dados;class(Y)
colnames(Y)<-c("North", "East", "South", "West") #Nomeando as colunas
Cork<- data.frame(North =y1, East =y2,South = y3, West =y4) #Criando um dataframe; class(Cork)


## Obtendo os resultados
p <- ncol(Y) #número de colunas
n <- nrow(Y) #número de linhas
In <- diag(n) # cria uma matriz identidade n x n
jn <- matrix(1, nrow = n, ncol = 1) # cria um vetor coluna de 1s com comprimento n
Jnn <- jn%*%t(jn) #matrix(1, nrow = n, ncol = n)# cria uma matriz n x n de 1s
Sigma <- (1 / (n - 1)) * t(Y) %*% (In - (1 / n) * Jnn) %*% Y #calcula a matriz de variâncias e covariâncias
D <- sqrt(diag(Sigma)) #Raiz quadrada da diagonal da matriz de variâncias e covariâncias, ou seja, raiz quadrada das variâncias
inv_D <- solve(diag(D)) # calcula a inversa da matriz diagonal D
corr <- inv_D %*% Sigma %*% inv_D # calcula a matriz de correlação
rownames(corr) <- c("North", "East", "South", "West") #Apenas nomeando as colunas
colnames(corr)<-c("North", "East", "South", "West")#Apenas nomeando as linhas
Verifica <- diag(sqrt(diag(Sigma)))%*% corr %*%diag(sqrt(diag(Sigma))) #Retorna a matriz de variâncias e covariâncias
rownames(Verifica) <- c("North", "East", "South", "West")#Apenas nomeando as linhas
colnames(Verifica)<-c("North", "East", "South", "West")#Apenas nomeando as colunas

##
#Imprimindo os resultados
cat("Matriz de variâncias e covariâncias amostrais:\n")
print(Sigma, digits = 4)

cat("\nMatriz de correlações amostrais:\n")
print(corr, digits = 5)

cat("\nVerificação: retorna a matriz de variâncias e covariâncias amostrais\n")
print(Verifica, digits = 4)
##

mi <- (1/n) * t(jn) %*% Y # Calcula o vetor de médias de cada coluna; colMeans(Y)
#Imprimindo os resultados
cat("Vetor de médias: ")
print(format(mi, digits=4), quote=FALSE) 



DM2 <- rep(0, n) # Inicia o vetor para o cálculo da distância de Mahalanobis
i <- 1 # Inicia o índice do loop while
while (i <= n) {
  yi <- Y[i, ] # Seleciona a linha i da matriz Y
  DM <- (yi - mi) %*% solve(Sigma) %*% t(yi - mi) # Calcula a distância de Mahalanobis
  DM2[i] <- DM # Atribui o valor calculado ao vetor de distâncias
  i <- i + 1 # Incrementa o índice
}


rank <- rank(DM2) # Calcula o rank das distâncias de Mahalanobis


## Apenas 'imprimindo' os resultados obtidos com a programação acima
# cat("-----------------------------------------------------------------\n",
#     "Distância de Mahalanobis de cada ponto (y) ao vetor de médias(mi)\n",
#     "-----------------------------------------------------------------\n")
# print(paste(Y, " ", format(DM2, digits=4), " ", rank, sep="")) # para cada valor em relação ao vetor de médias


cat("----------------------------------------------------------------------\n",
    "Distância de Mahalanobis de cada ponto (Y[i,]) ao vetor de médias(mi)\n",
    "---------------------------------------------------------------------\n")
print(cbind(Y, DM2, rank)) # considera a linha completa (vetores[1:28,]) em relação ao vetor de médias (mi)


#Usando funções implemenadas no R e imprimindo os resultados obtidos
cov_matrix <- cov(Cork[, c('North', 'East', 'South', 'West')]) # Calcula a matriz de variâncias e covariâncias com uma função implementada
cat("-----------------------------------------------------------------\n",
    'Matriz de variâncias e covariâncias utilizando a função cov():\n',
    "-----------------------------------------------------------------\n")
print(cov_matrix) # Imprime a matriz 

cor_matrix <- cor(Cork[, c('North', 'East', 'South', 'West')]) # Calcula a matriz de correlação com uma função implementada no R
cat("-----------------------------------------------------------------\n",
    'Matriz de correlação utilizando a função cor():\n',
    "-----------------------------------------------------------------\n")
print(cor_matrix) # Imprime a matriz 

d.m <- mahalanobis(Cork, colMeans(Cork), cov(Cork))
cat("-----------------------------------------------------------------------------------------------\n",
    'Distância de Mahalanobis de cada vetor (Y[i,]) ao vetor de médias usando a função mahalanobis()\n',
    "----------------------------------------------------------------------------------------------\n")
print(d.m) # 'Imprime' o vetor 


verifica.mahalanobis <- cbind(DM2,d.m)
colnames(verifica.mahalanobis) <- c("Função programada","Função implementada no R") #Verifica os valores da função programada e da função implementada no R

cat("----------------------------------------------------------------------------------------------------------\n",
    'Verificação dos valores de distância de Mahalanobis com a função programada e a função implementada no R\n',
    "---------------------------------------------------------------------------------------------------------\n")
print(verifica.mahalanobis)


# Aula 2 ------------------------------------------------------------------
#Construção do gráfico da normal bivariada no R 
rm(list = ls()) #limpando a memória do R
library(plotly) #carregando a biblioteca plotly


r <- 0 #Valor de correlação (altere esse valor [-1,1] e verifique o seu efeito no gráfico (3D dinâmico) da f.d.p. normal bivariada)
pi <- 3.1416 #Valor de pi (já é definido no R)
y1 <- seq(-4, 4, by = 0.1) #Definindo os valores da variável y1
y2 <- seq(-4, 4, by = 0.1) #Definindo os valores da variável y2
z <- outer(y1, y2, function(x,y) { 
  phi <- exp(-(x^2 - 2*r*x*y + y^2)/(2*(1-r^2))) / (2*pi*sqrt(1-r^2))
  return(phi)
}) #Função de densidade de probabilidade normal bivariada 

#Plotando o gráfico usando a biblioteca plotly
plot_ly(x = y1, y = y2, z = z, type = "surface") %>%
  layout(scene = list(xaxis = list(title = "y1"), 
                      yaxis = list(title = "y2"),
                      zaxis = list(title = "dens")),
         title = paste("Densidade Normal Bivariada (r =", r, ")"))



# Exemplo 4.4(a) ----------------------------------------------------------

rm(list = ls())
y <- c('y1','y2','y3')
mi <- c(3,1,2)
Sigma <- matrix(c(4, 0, 2, 0, 1, -1, 2, -1, 3), nrow = 3, ncol = 3)
##
a <- c(1,-2, 1)
mi_z <- t(a) %*% mi; mi_z #E(z)
var_z <- t(a) %*% Sigma %*% a; var_z #var(z)

##
ZZ <- matrix(c('z1 = y1-y2+y3','z2 = 3y1+y2-2y3'), nrow = 2, ncol = 1)
A <- matrix(c(1,-1,1, 3,1,-2),byrow = T, nrow = 2, ncol = 3)
mi_ZZ <- A %*% mi; mi_ZZ  ##E(A)
Sigma_ZZ <- A %*% Sigma %*% t(A); Sigma_ZZ #cov(Z)

##
A12 <- matrix(c(1,0,0, 0,1,0), byrow = T,nrow = 2, ncol = 3)
mi_12 <- A12 %*% mi;mi_12 #E(A12)
Sigma_12 <- A12 %*% Sigma %*% t(A12);Sigma_12 #cov(A12)

##
A13 <- matrix(c(1,0,0, 0,0,1), byrow = T,nrow = 2, ncol = 3)
mi_13 <- A13 %*% mi; mi_13 #E(A13)
Sigma_13 <- A13 %*% Sigma %*% t(A13); Sigma_13#cov(A13)

##
a1 <- as.matrix(c(1,0,0))
b2 <- as.matrix(c(0,1,0))
cov_12 <- t(a1) %*% Sigma %*% b2;cov_12  #cov(a1b2)


# Exemplo 4.4(b) ----------------------------------------------------------
rm(list = ls())
mi <- c(2,5,-2,1)
Sigma <- matrix(c(9,0,3,3,0,1,-1,2,3,-1,6,-3,3,2,-3,7),byrow = T, nrow=4)

##
Ay <- matrix(c(1,0,0,0,0,1,0,0),byrow = T, nrow=2)
(mi_y <- Ay %*% mi) #E(y)
(Sigma_yy <- Ay %*% Sigma %*% t(Ay)) #cov(y)

##
Ax <- matrix(c(0,0,1,0,0,0,0,1),byrow = T, nrow=2)
(mi_x <- Ax %*% mi)#E(x)
(Sigma_xx <- Ax %*%Sigma %*% t(Ax))#cov(x)
#
(Sigma_yx <- Ay %*% Sigma %*% t(Ax))#cov(yx)

(cov_ydx <- Sigma_yy - Sigma_yx %*% solve(Sigma_xx) %*% t(Sigma_yx))


#  Exemplo 4.4(c): -------------------------------------------------------
rm(list = ls())
#
mi <- c(2, 5, -2, 1)
Sigma <- matrix(c(9, 0, 3, 3, 0, 1, -1, 2, 3, -1, 6, -3, 3, 2, -3, 7), nrow=4, ncol=4, byrow=TRUE)

##
Ay <- matrix(c(1, 0, 0, 0), nrow=1, ncol=4)
(mi_y <- Ay %*% mi) # E(y)
(Sigma_yy <- Ay %*% Sigma %*% t(Ay))# cov(y)

##
Ax <- matrix(c(0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1), nrow=3, ncol=4, byrow=TRUE)
(mi_x <- Ax %*% mi) # E(x)
(Sigma_xx <- Ax %*% Sigma %*% t(Ax))# cov(x)

##
(Sigma_yx <- Ay %*% Sigma %*% t(Ax))# cov(yx)


(cov_ydx <- Sigma_yy - Sigma_yx %*% solve(Sigma_xx) %*% t(Sigma_yx))



# Exemplo 4.5 -------------------------------------------------------------
rm(list = ls())

Sigma <- matrix(c(9, 0, 3, 3, 0, 1, -1, 2, 3, -1, 6, -3, 3, 2, -3, 7), nrow = 4, byrow = TRUE)
D <- sqrt(diag(Sigma))
(Ro <- solve(diag(D)) %*% Sigma %*% solve(diag(D))) #cor matrix



Ay <- matrix(c(1, 0, 0, 0, 0, 1, 0, 0), byrow=T, nrow = 2)
(Sigma_yy <- Ay %*% Sigma %*% t(Ay))  #cov(y)

Dyy <- sqrt(diag(Sigma_yy))
(Ro_yy <- solve(diag(Dyy)) %*% Sigma_yy %*% solve(diag(Dyy)))#cor matrix(yy)

Ax <- matrix(c(0, 0, 1, 0, 0, 0, 0, 1), byrow=T, nrow = 2)
(Sigma_xx <- Ax %*% Sigma %*% t(Ax))#cov(x)

(Sigma_yx <- Ay %*% Sigma %*% t(Ax))#cov(yx)

(cov_ydx <- Sigma_yy - Sigma_yx %*% solve(Sigma_xx) %*% t(Sigma_yx))
D <- sqrt(diag(cov_ydx))
(Ro_ydx <- solve(diag(D)) %*% cov_ydx %*% solve(diag(D)))

Ro_yy
Ro_ydx



# Exemplo 6.2-CAPÍTULO 6. REGRESSÃO LINEAR SIMPLES ------------------------------------
###Exemplo 6.2

rm(list = ls())

cat("\nExemplo 6.2. Relation between exam score (y) and homework score (x)\n")

y <- as.vector(c(95,80,0,0,79,77,72,66,98,90,0,95,35,50,72,55,75,66))
x1 <- as.vector(c(96,77,0,0,78,64,89,47,90,93,18,86,0,30,59,77,74,67))
n <- length(y)
k <- 1
jn <- matrix(rep(1,n), ncol=1)
Jnn <- jn%*%t(jn)
In <- diag(n)
X <- matrix(cbind(jn, x1),ncol = 2)
y_barra <- (1/n)*Jnn%*%y # (1/n)*Jnn=MG
Tot <- y-y_barra
SQTotal <- t(Tot)%*%Tot

Beta1 <- as.numeric(t(x1)%*%(In-(1/n)*Jnn)%*%y/(t(x1)%*%(In-(1/n)*Jnn)%*%x1))

Beta0 <- as.numeric((1/n)*t(jn)%*%(y - Beta1*x1))

cat("\nEstimativas dos parâmetros da reta:\n")
cat("Beta0: ", Beta0, "\n")
cat("Beta1: ", Beta1, "\n")

lm(y~x1) #Função implementada no R
anova(lm(y~x1))

##
y_hat <- Beta0 + Beta1*x1
Reg <- y_hat - y_barra
SQReg <- t(Reg)%*%Reg
Res <- y - y_hat
SQRes <- t(Res)%*%Res
s2 <- as.numeric(SQRes/(n - k - 1))
# s2 <- (t(y)*y - t(Beta)*t(X)*y)/(n-k-1)
res_pad <- Res/sqrt(s2)



cat("Valores observados(y) e estimados(y_hat), residuo(res) e residuo padronizado(res_pad):\n",
    "--------------------------------------------------------------------------------------\n")
print(data.frame(y, y_hat, Res, res_pad), digits=4)

var_y <- (t(y) %*% (In - (1/n)*Jnn) %*% y)/(n-1) # Calcula a variância amostral de y
s <- sqrt(s2)

cat("Variância dos dados originais:", var_y, "\n",
    "Variância de y|x: ", s2, "\n",
    "Desvio padrão de y|x : ", s, "\n\n")

x_barra <- t(jn) %*% (x1/n);x_barra 
var_Beta0 <- s2*((1/n) + (x_barra^2/(t(x1) %*% (In - (1/n)*Jnn) %*% x1)));var_Beta0
stderr_Beta0 <- sqrt(var_Beta0)
var_Beta1 <- s2/(t(x1) %*% (In - (1/n)*Jnn) %*% x1);var_Beta1
stderr_Beta1 <- sqrt(var_Beta1)
vcov(lm(y~x1))#Função implementada no R

ttab <- qt(0.975, n-2) #Obtendo o quantil t
liminf0 <- Beta0 - ttab*stderr_Beta0; limsup0 <- Beta0 + ttab*stderr_Beta0
liminf1 <- Beta1 - ttab*stderr_Beta1; limsup1 <- Beta1 + ttab*stderr_Beta1

cat("Beta0\t Var(Beta0)\t Std.Err.(Beta0)\t IC 95% (Beta0)\n",
    Beta0, "\t", var_Beta0, "\t", stderr_Beta0, "\t[", liminf0, ", ", limsup0, "]\n",
    "Beta1\t Var(Beta1)\t Std.Err.(Beta1)\t IC 95% (Beta1)\n",
    Beta1, "\t", var_Beta1, "\t", stderr_Beta1, "\t[", liminf1, ", ", limsup1, "]\n\n")

R2 <- SQReg/SQTotal # Coeficiente de determinação - R2
corr <- sqrt(R2)

cat("SQTotal = SQReg + SQRes\n",
    SQTotal, "\t", SQReg, "\t", SQRes, "\n\n",
    "Coeficiente de determinação (R2): ", R2, "\n",
    "Coeficiente de correlação (r): ", corr, "\n\n")

(tcalc1 <- Beta1/stderr_Beta1) # Para testar H0: Beta1 = 0
(tcalc2 <- corr*sqrt(n-2)/(sqrt(1-corr^2)))
(p_valor <- 2*(1-pt(abs(tcalc1),n-2)))

# cat("H0: Beta1 = 0\t", format(tcalc1, digits = 4), 
#     format(tcalc2, digits = 4), format(p_valor, digits = 4), "\n")


# Definição dos dados
dados <- data.frame(Ind = 1:18, 
                    y_exam = c(95, 80, 0, 0, 79, 77, 72, 66, 98, 90, 0, 95, 35, 50, 72, 55, 75, 66),
                    x_homework = c(96, 77, 0, 0, 78, 64, 89, 47, 90, 93, 18, 86, 0, 30, 59, 77, 74, 67))

# Análise de regressão
modelo <- lm(y_exam ~ x_homework, data = dados)

# Imprime resultados do teste de hipótese H0: Beta1 = 0
cat("H0: Beta1 = 0\t")
(tcalc1 <- coef(modelo)[2] / summary(modelo)$coef[2, 2])
(tcalc2 <- qt(0.025, df = summary(modelo)$df))
(p_valor <- 2 * pt(abs(tcalc1), df = summary(modelo)$df, lower.tail = FALSE))

# cat("H0: Beta1 = 0\t", format(tcalc1, digits = 4), 
#     format(tcalc2, digits = 4), format(p_valor, digits = 4), "\n")

# Imprime a saída da análise de regressão
print(summary(modelo))

# Plot dos dados e ajuste da reta de regressão
library(ggplot2)
ggplot(dados, aes(x = x_homework, y = y_exam)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE) +
  labs(x = "x_homework", y = "y_exam") 


# REGRESSÃO LINEAR MÚLTIPLA: ESTIMAÇÃO  -----------------------------------
#Exemplo 7.2

#library(Matrix)
rm(list = ls())
y <- c(2,3,2,7,6,8,10,7,8,12,11,14)
n <- length(y)
In <- diag(n)
jn <- matrix(rep(1,n),ncol = 1)
Jnn <-jn %*% t(jn)

x0 <- jn
x1 <- c(0,2,2,2,4,4,4,6,6,6,8,8)
x2 <- c(2,6,7,5,9,8,7,10,11,9,15,13)

k <- 2

X01 <- cbind(x0, x1)
X02 <- cbind(x0, x2)
X012 <- cbind(x0, x1, x2)

Beta01 <- solve(t(X01) %*% X01) %*% t(X01) %*% y
Beta02 <- solve(t(X02) %*% X02) %*% t(X02) %*% y
Beta012 <- solve(t(X012) %*% X012) %*% t(X012) %*% y

cat("  Modelo: y = b0 + b1*x1 + e          =>     Beta = ", format(Beta01, digits=4), "\n")
cat("  Modelo: y = b0 + b2*x2 + e          =>     Beta = ", format(Beta02, digits=4), "\n")
cat("  Modelo: y = b0 + b1*x1 + b2*x2 + e  =>     Beta = ", format(Beta012, digits=4), "\n")

Beta01
Beta02
Beta012

#cat(paste(format(Beta01, digits=4), format(Beta02, digits=4), format(Beta012, digits=4)))


X <- cbind(x0, x1, x2);colnames(X) <-c("x0", "x1", "x2") 
Beta <- round(solve(t(X) %*% X) %*% t(X) %*% y,4)
y_hat <- round(X %*% Beta,4); colnames(y_hat) <- ("y_hat")
Beta
(XXX <- cbind(X,y,y_hat))


####página 168 (Modelo centrado)
x1x2 <- cbind(x1, x2)
x1x2c <- round((In - (1/n)*Jnn)%*% x1x2,4) #obtendo os valores de x's centrados
Xc <- cbind(x0, x1x2c);colnames(Xc) <-c("x0", "x1c", "x2c")
Betac <- round(solve(t(Xc) %*% Xc) %*% t(Xc) %*% y,4)#obtendo as estimativas dos coeficientes de regressão com x's centrados
y_hatc <- round(Xc %*% Betac,4);colnames(y_hatc) <- ("y_hatc")

Betac
XXXc <-cbind(Xc,y,y_hatc)

xxunido <- list(Beta,XXX,Betac,XXXc); names(xxunido) <- c("Beta","XXX","Betac","XXXc");xxunido #unindo os resultados



XcLXc <- t(Xc) %*% Xc;XcLXc


###Pagina 175
res <- y - y_hatc
SQRes <- t(res) %*% res
s2 <- as.numeric((1/(n-k-1)) * SQRes)
cov_Beta <- round(s2 * solve(t(Xc) %*% Xc),4);
data.joing <- list(Betac,s2,cov_Beta);names(data.joing) <-c("Betac","s2","cov_Beta");data.joing

####Pagina 195
(SQreg <- t(y) %*% (X %*% solve(t(X) %*% X) %*% t(X) - (1/n) * Jnn) %*% y)
(SQTot <- t(y) %*% (In - (1/n) * Jnn) %*% y)
(R2 <- SQreg / SQTot)
(R2aj <- ((n-1) * R2 - k) / (n - k - 1))


# Exemplo 7.2 - Ortogonalizacao -------------------------------------------

rm(list = ls())
y <- c(2, 3, 2, 7, 6, 8, 10, 7, 8, 12, 11, 14)
n <- length(y)
jn <- matrix(1, n, ncol=1)
x1 <- c(0, 2, 2, 2, 4, 4, 4, 6, 6, 6, 8, 8)
x2 <- c(2, 6, 7, 5, 9, 8, 7, 10, 11, 9, 15, 13)
X <- cbind(jn, x1, x2);colnames(X) <- c("b0", "x1", "x2")
X01 <- cbind(jn, x1);colnames(X01) <- c("b0", "x1")

#1) Ajusta a regressao de y sobre x1;
Beta_yx0x1 <- solve(t(X01) %*% X01) %*% t(X01) %*% y #lm(y~x1)
cat("Beta_yx0x1: ", format(Beta_yx0x1, digits = 4), "\n")
y_hatyx0x1 <- X01 %*% Beta_yx0x1 #lm(y~x1)$fitted.values
Res_yx0x1 <- y - y_hatyx0x1 #Calcula o residuo da regressao de y sobre x1;
org.data <- data.frame(y,y_hatyx0x1,Res_yx0x1)
xx <- list(Beta_yx0x1,org.data);names(xx) <- c("Beta_yx0x1","org.data");xx



#2) Ajusta a regressao de x2 sobre x1
Beta_x2x1 <- solve(t(X01) %*% X01) %*% t(X01) %*% x2
Beta_x2x1
x2_hatx2x1 <- X01 %*% Beta_x2x1
x2_hatx2x1
Res_x2x1 <- x2 - x2_hatx2x1 #Calcula o residuo da regressao de x2 sobre x1
Res_x2x1
org.data2 <- data.frame(x2,x2_hatx2x1,Res_x2x1)
xx2 <- list(Beta_x2x1,org.data2);names(xx2) <- c("Beta_yx2x1","org.data2");xx2



#3) Ajusta a regressao de Res_yx0x1 sobre Res_x2x1;

Beta2 <- solve(t(Res_x2x1) %*% Res_x2x1) %*% t(Res_x2x1) %*% Res_yx0x1
Beta2 

table <- cbind(y, x1, x2, round(Res_yx0x1,4),  round(Res_x2x1,4))
colnames(table) <- c("y", "x1", "x2", "Res_yx0x1", "Res_x2x1")
table #Tabela 7.2 - Dados da Tabela 7.1 e residuos


Beta <- round(solve(t(X) %*% X) %*% t(X) %*% y,4)
Beta


verifica <- t(Res_x2x1) %*% X01
verifica #Verifica a ortogonalidade de Res_x2x1 e as colunas de X01



# Exemplo 7.8.2(a) --------------------------------------------------------


#library(Matrix) # Carrega a biblioteca Matrix para usar a função 'sparseMatrix()'
rm(list = ls())
y <- c(2,5,5,8,6,9,7,11,7,15,10,17)
x1 <- c(2,5,6,7,7,8,9,9,10,11,13,15)
n <- length(y)
x0 <- matrix(1, nrow = n, ncol = 1)


X <- cbind(x0, x1)
k <- ncol(X) - 1
Beta <- solve(t(X) %*% X) %*% t(X) %*% y # Aqui usamos a função solve() para obter a matriz inversa de (t(X) %*% X)
y_hat <- X %*% Beta
res <- y - y_hat
s2 <- as.numeric((1/(n-k-1)) * t(res) %*% res)
cov_Beta <- s2 * solve(t(X) %*% X)
cat("Modelo com estrutura cov(y) = I*sigma2\n\n")
Beta 
s2 
cov_Beta


V <- diag(x1)
Betaz <- round(solve(t(X) %*% solve(V) %*% X) %*% t(X) %*% solve(V) %*% y,4)
y_hatz <- round(X %*% Betaz,4)
resz <- y - y_hatz
s2z <- as.numeric((1 / (n - k - 1)) * t(resz) %*% solve(V) %*% resz)
cov_Betaz <- s2z * solve(t(X) %*% solve(V) %*% X)
print("Modelo com estrutura cov(y) = V*sigma2")
Betaz
s2z
cov_Betaz



##
library(data.table)
K <- cbind(x1, y, y_hat, y_hatz)
colnames(K) <- c("x", "y", "y_hat", "y_hatz")
graf <- data.table(K)


library(ggplot2)
ggplot(graf, aes(x = x, y = y)) +
  geom_point() +
  geom_line(aes(y = y_hat), color = "blue",show.legend = "y_hat") +
  geom_line(aes(y = y_hatz), color = "red",show.legend = "y_hatz") +
  ggtitle("Ajustes MQ Ordinario (y_hat) e Ponderado (y_hatz)")