---
title: "New orders in manufacturing sector"
author: "Hele-Liis Peedosk"
date: "18 aprill 2019"
output: html_document
---

# Load libraries and functions. Run h2o.
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r libraries} 
source("functions.R")
load_libraries()
```

```{r}
# Run 'h2o' interface
init_h2o(install=FALSE)
# h2o.shutdown(prompt=FALSE)
```


# Initialize results table
```{r}
results=init_results()
```

## Read in the data

```{r}    
data_dir <- "C:/Users/helts/Desktop/Masters thesis/R"
fname <- file.path(data_dir, "USM_monthly.csv")
data_USM = read.csv(fname, header=TRUE, sep=";", dec=",",na.strings="0") 
data_USM = data.matrix(data_USM[,-1], rownames.force = NA)
data_USM=as.data.frame(matrix(c(t(data_USM)), ncol=1))
colnames(data_USM)=c("value")
data_USM$date=seq(as.Date("1958/1/1"), by = "month", length.out = 12*(2019-1958))

# Change the date from factor type to date type, add columns month, weekday.
dates<-data_USM$date
```

```{r}
# Create tibble table with differenced values
trend=FALSE
diff=TRUE
data=data_USM[data_USM$date>'1979-12-31',]
tbl_data=create_tbl(data_USM[data_USM$date>'1979-12-31',], remove_trend=trend, difference=diff)
USM=tbl_data$data
data=tbl_data$tbl
data_orig=tbl_data$tbl_orig
data_trend=tbl_data$tbl_trend

# To have the first value for the differenced data
data$value[1]=145099-141003
```

# Draw a graph of the time series
```{r}
source("functions.R")

df_orig=train_val_test_split(data_orig,val_split=val_split, test_split = 0.85)
latex_graph_df(series=df_orig, subseries=df_orig[df_orig$index>'2009-12-31',], y_lab="Millions of dollars", y_axis_lab_scale=scales::comma)
```


# Define the vaiables for this time series
```{r}
seasonality=12

# Choose input size, horizon, hidden layers, seed, epochs, activation_function
seed = 1
epochs = 50
activation_function = "Rectifier"

# Information for results table
data_name = "USM"
rolling = FALSE
comments = "-"
history_start = "Starting from 1980"


# Choose whether to use validation set
val_split = 0.8
val_set = TRUE

# ARIMA orders
order=c(0,1,1)
seasonal=list(order=c(0,1,1), period=seasonality)

```



# Determine manual grid

```{r}
# Grid search parameters
inputs= c(12, 24, 36)
horizons = c(1, 3, 6, 12)        #c(1, 3, 5, 10, 20, 60)   # daily
hidden_layers = list(c(3),c(5), c(8), c(10), c(10, 10), c(20, 20), c(100),c(200), c(300))
ensemble_sizes = c(10) 

```


# Develop neural network models 

```{r}

df=train_val_test_split(data,val_split=val_split, test_split = 0.85)
df_orig=train_val_test_split(data_orig,val_split=val_split, test_split = 0.85)

for (n_output in horizons){
    for (n_input in inputs){    
       tryCatch({      # In case the input horizon is too long.
         
          # Generate data with certain colnames, extract data, X colnames, y colnames
          io_df<-gen_input_output(df, n_input, n_output, val_set=TRUE, step_fn=NULL, key_val=NULL)
          io_df_orig<-gen_input_output(df_orig, n_input, n_output, val_set=TRUE, step_fn=NULL, key_val=NULL, standardize=FALSE)

          for (ensemble_size in ensemble_sizes){

              for (hidden in hidden_layers){
                  if (sum(hidden)>=20) model_name="h2o: large" else model_name="h2o: small"
                  tryCatch({      # In case of some error with the model
                  train_eval(io_df=io_df, ensemble_size=ensemble_size, hidden=hidden,  data_name=data_name, train_predict_model = train_predict_h2o_sep
                             , package="h2o", comments="multi NN", rolling=rolling, history_start=history_start
                             , mod_fn=mod_fn_h2o_reg, model_name=model_name, print_info=TRUE, hyper_params=NULL, io_df_orig=io_df_orig
                             , write_results_path= "", diff=diff, trend=trend)
                  },error=function(e){cat("ERROR: Cannot handle the h2o model with horizon=",n_output
                                       , ", lookback=",n_input, ", and hidden ", paste0(hidden, collapse = ", "), "\n")})


                  tryCatch({      # In case of some error with the model
                  train_eval(io_df=io_df, ensemble_size=ensemble_size, hidden=hidden,  data_name=data_name, train_predict_model = train_predict_h2o_rec
                             , package="h2o", comments=paste(comments, "iterative", sep=", "), rolling=rolling, history_start=history_start
                             , mod_fn=mod_fn_h2o_reg, model_name=model_name, print_info=TRUE, hyper_params=NULL, io_df_orig=io_df_orig
                             , write_results_path= "", diff=diff, trend=trend)
                  },error=function(e){cat("ERROR: Cannot handle the h2o model with horizon=",n_output
                                       , ", lookback=",n_input, ", and hidden ", paste0(hidden, collapse = ", "), "\n")})

                  if (sum(hidden)>=20)  tryCatch({      # In case of some error with the model
                      train_eval(io_df=io_df, ensemble_size=ensemble_size, hidden=hidden,  data_name=data_name, train_predict_model = train_predict_neuralnet
                             , package="neuralnet", comments="direct", rolling=rolling, history_start=history_start, diff=diff, trend=trend
                             , mod_fn=mod_fn_h2o_reg, model_name="small", print_info=TRUE, hyper_params=NULL, io_df_orig=io_df_orig
                             , write_results_path= "")
                  },error=function(e){cat("ERROR: Cannot handle the neuralnet model with horizon=",n_output
                                           , ", lookback=",n_input, ", and hidden ", paste0(hidden, collapse = ", "), "\n")})
 


               }}

          tryCatch({      # In case of some error with the model
                  train_eval(io_df=io_df_orig, 1, c(1),  data_name, train_predict_arima, package="ARIMA", comments, rolling, history_start
                         , mod_fn=NULL, model_name="SARIMA", print_info=TRUE, hyper_params=NULL,  seasonality=seasonality, io_df_orig=io_df_orig
                         , order=order, seasonal=seasonal, frequency=frequency
                         , diff=diff, trend=trend, val_set=TRUE
                         , write_results_path= "")


          
          tryCatch({      # In case of some error with the model
                  train_eval(io_df=io_df_orig, ensemble_size, hidden,  data_name, train_predict_rwf, package="Naive", comments, rolling, history_start
                         , mod_fn=NULL, model_name="snaive", print_info=TRUE, hyper_params=NULL, seasonality=seasonality, io_df_orig=io_df_orig
                         , write_results_path=data_dir)
          },error=function(e){cat("ERROR: Cannot handle naive model with horizon =",n_output
                               , "and lookback=",n_input, "\n")})



          },error=function(e){cat("ERROR: Cannot handle arima model with horizon =",n_output
                               , "and lookback=",n_input, "\n")})
          # 
      },error=function(e){cat("ERROR: Cannot handle input-ouput split with horizon =",n_output
                           , "and lookback=",n_input, "\n")})

  }}

    results
```



# ARIMA for input size 12
We need to take additional data from train/validation sets, as SARIMA needs more input as 12. 
```{r}
plot(df_orig$value)
```


```{r}
inputs= c(12)
horizons = c(1, 3, 6, 12)        #c(1, 3, 5, 10, 20, 60)   # daily
fact=2
order=c(1,0,0)
seasonal=list(order=c(0,1,2), period=12)

df=train_val_test_split(data,val_split=val_split, test_split = 0.85)
df_orig=train_val_test_split(data_orig,val_split=val_split, test_split = 0.85)


for (n_output in horizons){
      for (n_input in inputs){    
       tryCatch({      # In case the input horizon is too long.
         
          # Generate data with certain colnames, extract data, X colnames, y colnames
          io_df<-gen_input_output(df, n_input, n_output, val_set=TRUE, step_fn=NULL, key_val=NULL)
          io_df_orig<-gen_input_output(df_orig, n_input, n_output, val_set=TRUE, step_fn=NULL, key_val=NULL, standardize=FALSE)        

          df_valid=df_orig
          df_test=df_orig

          df_valid[(sum(df_valid$key=="training")-seasonality*fact+1):sum(df_valid$key=="training"),"key"]="validation"
          df_test[(sum(df_test$key=="validation")+sum(df_test$key=="training")-seasonality*fact+1):(sum(df_test$key=="validation")
                                                                                                    +sum(df_test$key=="training")),"key"]="testing"


          io_df_valid<-gen_input_output(df_valid, n_input+seasonality*fact, n_output, val_set=TRUE, step_fn=NULL, key_val="validation"
                                        , standardize=FALSE, center_history = 0
                                        , scale_history = 1)
          io_df_test<-gen_input_output(df_test, n_input+seasonality*fact, n_output, val_set=TRUE, step_fn=NULL, key_val="testing"
                                        , standardize=FALSE, center_history = 0
                                        , scale_history = 1)
          
          io_df_orig$X_valid=io_df_valid$X_valid
          io_df_orig$X_test=io_df_test$X_test

          
          tryCatch({      # In case of some error with the model
          train_eval(io_df=io_df_orig, 1, c(1),  data_name, train_predict_arima, package="ARIMA", comments, rolling, history_start
                         , mod_fn=NULL, model_name="SARIMA", print_info=TRUE, hyper_params=NULL, seasonality=seasonality, io_df_orig=io_df_orig
                         , order=order, seasonal=seasonal
                         , diff=TRUE, trend=FALSE, val_set=TRUE
                         , write_results_path= data_dir)

          },error=function(e){cat("ERROR: Cannot handle arima model with horizon =",n_output
                               , "and lookback=",n_input, "\n")})
          
      },error=function(e){cat("ERROR: Cannot handle input-output split with horizon =",n_output
                           , "and lookback=",n_input, "\n")})

  }}



    results
```



# Autocorrelation and partial autocorrelation plots for the series


```{r}
Z1<-ts(df_orig$value[df$key=="training"],frequency = 12, start=decimal_date(df_orig$index[1]))
plot(Z1)

par(family = "LM Roman 10", mfrow=c(2,1))
acf(Z1,32, main = "")
title(sub="(a) Autocorrelations of electricity market prices.", adj = 0)
pacf(Z1,32, main="")
title(sub="(b) Partial autocorrelations of electricity market prices.", adj = 0)
```


