install.packages(c("ARDL", "lmtest", "normtest"))
library(ARDL)
#library(urca)
library(lmtest)
library(normtest)

Sys.setenv(`_R_S3_METHOD_REGISTRATION_NOTE_OVERWRITES_` = "false")


options(stringsAsFactors = F)
# I combined Gretl with R and Microfit

# For the ADF unit root test I combined R and Gretl

# Loading the database to obtain the ARDL estimations
cuba_ardl <- read.csv("~/Desktop/Website/gdp_infra.csv")
cuba_ardl <- ts(cuba_ardl[-1], start = 1920, frequency = 1)

# in the file all outputs.txt are all the results of the ARDL systems
#sink(file = "all outputs.txt")

######## Best ARDL model selection #####

AIC_selection <- auto_ardl(data = cuba_ardl, max_order = 5, selection = "AIC",  
                           formula = gdp ~ energy + train | crisis29, 
                           selection_minmax = "min",  search_type = "horizontal",
                           start = 1920, end = 1957, grid = TRUE)
#   fixed_order = -1, starting_order = NULL, 
AIC_selection

######## Model estimation #####

ardl_model<- ardl(data = cuba_ardl, order = AIC_selection$best_order, 
                  formula = gdp ~ energy + train | crisis29, start = 1920, end = 1957)
summary(ardl_model)

######## Breusch-Godfrey test upto 4th order #####

for (i in 1:2){
    godfrey<-bgtest(ardl_model, order = i, type = "Chisq")
    print(godfrey)
}

######## Breusch-Pagan test #####
print(bptest(ardl_model, studentize = TRUE))

######## Shapiro Wilk test #######

shapiro.test(ardl_model$residuals)

######## Jarque Bera test #######

jb.norm.test(ardl_model$residuals)

######## adjusted Jarque Bera test (Urzua, 1996) #######

ajb.norm.test(ardl_model$residuals)

##### Residual plotting

plot(ardl_model$residuals,)

######## Ramsey's RESET test ######

for (i in 2:3){
    resettestt<-resettest(ardl_model, power = 2:i, type = "fitted")
    print(resettestt)
}

######## Case ######

print("case III")

######## F Bounds test ######

f_bound <- bounds_f_test(ardl_model, case = 3, exact = TRUE, R = 40000, test = "F" #, vcov_matrix = NULL
)
f_bound

######## t Bounds test ######

t_bound <- bounds_t_test(ardl_model, case = 3, exact = TRUE, R = 40000)
t_bound

######## Long-run equation ######
long_run <- multipliers(ardl_model, type = "lr", vcov_matrix = NULL)
print(long_run)

######## Unrestricted ECM equation ######
restricted<- recm(ardl_model, case = 3)
summary(restricted)