data('RecreationDemand')
h <- RecreationDemand
x <- h$trips
y <- h$income

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 tra khuyết tật bằng phương pháp Breusch-Pagan
  bp_test <- bptest(model)
  p_value <- bp_test$p.value
  
  # 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, p_value = p_value, 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 
## -2.8925 -0.8925 -0.4508  1.1075  5.1959 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.89246    0.07652  50.869   <2e-16 ***
## x           -0.01767    0.01146  -1.541    0.124    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.85 on 657 degrees of freedom
## Multiple R-squared:  0.003604,   Adjusted R-squared:  0.002087 
## F-statistic: 2.376 on 1 and 657 DF,  p-value: 0.1237
# Xem kết quả kiểm tra khuyết tật
result$p_value
##         BP 
## 0.01112883
# Xem kết quả dự báo
head(result$prediction, 10)
##        1        2        3        4        5        6        7        8 
## 3.892459 3.892459 3.892459 3.892459 3.892459 3.892459 3.892459 3.892459 
##        9       10 
## 3.892459 3.892459

HHHHHHHHHHHHHHHHH