library(AER)
data("Affairs")
h <- Affairs
x <- h$rating
y <- h$age
z <- h$occupation
regression_heteroskedasticity_forecast <- function(data, y_var, x_var) {
  # Chạy hồi quy tuyến tính
  model <- lm(formula = paste(y_var, "~", x_var), data = data)

  # Kiểm định tương quan
  cor <- cor.test(x,y)
  # Vẽ đồ thị
  plot <- plot(model)
  # Dự báo giá trị của biến phụ thuộc
  prediction <- predict(model, newdata = data)

  # Trả về kết quả
  return(list(model = model, cor = cor,plot = plot, prediction = prediction))
}
# Sử dụng function
result <- regression_heteroskedasticity_forecast(data.frame(x, y), "y", "x")

# Xem kết quả hồi quy
summary(result$model)
## 
## Call:
## lm(formula = paste(y_var, "~", x_var), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.549  -7.049  -2.049   4.627  26.302 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39.0755     1.3767   28.38  < 2e-16 ***
## x            -1.6756     0.3372   -4.97 8.76e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.111 on 599 degrees of freedom
## Multiple R-squared:  0.0396, Adjusted R-squared:  0.038 
## F-statistic:  24.7 on 1 and 599 DF,  p-value: 8.759e-07
# Xem kết quả kết quả độ tin cậy
result$cor
## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = -4.9698, df = 599, p-value = 8.759e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2746071 -0.1209471
## sample estimates:
##        cor 
## -0.1989999
# Xem kết quả dự báo của 6 giá trị đầu tiên
head(result$prediction)
##        1        2        3        4        5        6 
## 32.37321 32.37321 32.37321 30.69764 34.04879 30.69764