Market Segments

For each hotel, we have several market segments as mentioned earlier in the data description. In order to maximize our revenue, we will forecast the most profitable market segment. First, we will look at the number of transactions for each market segment on each hotel.

hotel %>% 
  count(hotel, market_segment, is_canceled) %>% 
  group_by(hotel) %>% 
  mutate(total = sum(n),
         ratio = n/total,
         is_canceled = ifelse(is_canceled == 0, "No", "Yes")) %>% 
  ungroup() %>% 
  mutate(market_segment = tidytext::reorder_within(market_segment, ratio, hotel)) %>%  
  ggplot(aes(ratio, market_segment, fill = is_canceled)) +
  geom_col(position = "dodge") +
  labs(x = "Percentage", y = "Market Segment", 
       title = "Hotel Demand by Market Segment",
       fill = "Booking Canceled") +
  facet_wrap(~hotel, scales = "free_y") +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  scale_fill_manual(values = c("skyblue", "firebrick")) +
  tidytext::scale_y_reordered() +
  theme_pander() +
  theme(legend.position = "top")

Based on the result, most of the booking done via travel agent, either online or offline, which combined together contributes more than 40% of the non-canceled total transactions. The other segment don’t have much transactions, but perhaps we would like to see the revenue generated by each market segment by looking at the Average Daily Rate (ADR). Average Daily Rate as defined by dividing the sum of all lodging transactions by the total number of staying nights. Therefore, the higher ADR means more revenue generated for each staying night.

hotel %>% 
  filter(is_canceled == F) %>% 
  group_by(hotel, market_segment) %>% 
  summarise(adr = sum(adr)) %>% 
  mutate(total = sum(adr),
         ratio = adr/total) %>% 
  ungroup() %>% 
  mutate(market_segment = tidytext::reorder_within(market_segment, ratio, hotel)) %>% 
  ggplot(aes(ratio, market_segment)) +
  geom_col(fill = "firebrick") +
  facet_wrap(~hotel, scales = "free_y") +
  tidytext::scale_y_reordered() +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  theme_pander() +
  labs(x = NULL, y = NULL, 
       title = "ADR Contributions by Market Segments")

Total ADR generated from online travel agent is the highest both in the city hotel and the resort hotel. The segment of offline travel agent and direct has a little margin with each other, with direct segment has the higher contribution in resort hotel even though it has lower number of transactions. This give us the top 3 market segments, both in term of quantity and profitability. We will focus on this segments for the rest of the analysis.

Time Series Analysis

We want to forecast the demand of lodging for both city hotel and resort hotel. Based on the previous exploration, we will use only the data from segment online travel agent, offline travel agent and direct. We will consider both canceled and non-canceled transactions to reflect the demand.

hotel_city <- hotel %>%
  filter(hotel == "City Hotel",
         market_segment %>% str_detect("TA|Direct")) 

hotel_resort <- hotel %>%
  filter(hotel == "Resort Hotel",
         market_segment %>% str_detect("TA|Direct")) 

hotel_agg_city <- hotel_city %>% 
  group_by(arrival_date) %>% 
  summarise(demand = n())

hotel_agg_resort <- hotel_resort %>% 
  group_by(arrival_date) %>% 
  summarise(demand = n())

head(hotel_agg_city, 10)

#> # A tibble: 10 × 2
#>    arrival_date demand
#>    <date>        <int>
#>  1 2015-07-01       79
#>  2 2015-07-02        8
#>  3 2015-07-03        3
#>  4 2015-07-04        6
#>  5 2015-07-05        8
#>  6 2015-07-06       10
#>  7 2015-07-07       18
#>  8 2015-07-08        8
#>  9 2015-07-09        6
#> 10 2015-07-10        7

Let’s look at the series of demand for each hotel.

hotel_agg_city   %>%
  mutate(hotel = "City Hotel") %>% 
  bind_rows(hotel_agg_resort %>% mutate(hotel = "Resort Hotel")) %>% 
  ggplot(aes(arrival_date, demand)) +
  geom_line(color = "dodgerblue4") +
  theme_pander() +
  facet_wrap(~hotel, ncol = 1) +
  labs(x = NULL, y = "Demand", title = "Hotel Booking Demand")

The demand for city hotel has a higher fluctuation compared to the resort hotel. This may be caused by several factors, including the room capacity, since we don’t know the room capacity for each hotel. There are also some spikes in demand for city hotel during the late 2015. We will inspect further by dividing the series by market segment. We will also need to make the series have constant interval of time, which is a daily interval in our case.

# summarizing data
hotel_agg_city <- hotel_city %>% 
  group_by(arrival_date, market_segment) %>% 
  summarise(demand = n())

hotel_agg_resort <- hotel_resort %>% 
  group_by(arrival_date, market_segment) %>% 
  summarise(demand = n())

# join city hotel and resort hotel
hotel_agg <- hotel_agg_city %>%
  mutate(hotel = "City Hotel") %>% 
  bind_rows(hotel_agg_resort %>% mutate(hotel = "Resort Hotel")) %>% 
  ungroup()

# padding
start_interval <-  ymd(range(hotel$arrival_date)[1])
end_interval <- ymd(range(hotel$arrival_date)[2])

hotel_pad <- hotel_agg %>% 
  group_by(hotel, market_segment) %>% 
  pad(start_val = start_interval, end_val = end_interval) %>% 
  replace_na(list(demand = 0)) %>% 
  ungroup()

hotel_pad %>% 
  ggplot(aes(arrival_date, demand, color = market_segment)) +
  geom_line() +
  theme_pander() +
  theme(legend.position = "top") +
  facet_wrap(~hotel, ncol = 1, scales = "free_y") +
  labs(x = NULL, y = "Demand", title = "Hotel Demand", color = NULL) +
  scale_color_manual(values = c("firebrick", "skyblue", "orange"))

The behaviour for each segment is quite different, so we will forecast them separately. At this point, we will have 6 different time series to be forecasted for 3 different segments on each hotel.

Seasonality Analysis

Before we proceed to create a forecasting model. We will try to gain more insight regarding the customer behaviour by looking at the seasonality of the demand. Since we have too many series, we will only explore the most lucrative market segment, that is the segment of online TA. We will look at the behaviour for both the city hotel and the resort hotel.

city_weekly <- hotel_agg_city %>% 
  filter(market_segment == "Online TA") %>% 
  ts(frequency = 7) %>% 
  decompose()

city_weekly$seasonal %>% 
  matrix(ncol = 7, byrow = T) %>%
  t() %>% 
  as.data.frame() %>% 
  select(V1) %>% 
  mutate(day = wday(head(hotel_agg_city %>% filter(market_segment == "Online TA") %>% pull(arrival_date) , 7), label = T)) %>% 
  ggplot(aes(day, V1, fill = V1)) +
  geom_col(show.legend = F) +
  geom_hline(yintercept = 0 ) +
  labs(x = NULL, y = "seasonality", title = "Weekly Seasonality of City Hotel, Online TA Segment ") +
  scale_fill_gradient2(low = "firebrick4", mid = "skyblue", high = "dodgerblue4") +
  theme_pander()

resort_weekly <- hotel_agg_resort %>% 
  filter(market_segment == "Online TA") %>% 
  ts(frequency = 7) %>% 
  decompose()

resort_weekly$seasonal %>% 
  matrix(ncol = 7, byrow = T) %>%
  t() %>% 
  as.data.frame() %>% 
  select(V1) %>% 
  mutate(day = wday(head(hotel_agg_resort %>% filter(market_segment == "Online TA") %>% pull(arrival_date) , 7), label = T)) %>% 
  ggplot(aes(day, V1, fill = V1)) +
  geom_col(show.legend = F) +
  geom_hline(yintercept = 0 ) +
  labs(x = NULL, y = "seasonality", title = "Weekly Seasonality of Resort Hotel, Online TA Segment ") +
  scale_fill_gradient2(low = "firebrick4", mid = "skyblue", high = "dodgerblue4") +
  theme_pander()

Monthly Seasonality

We will also check the monthly seasonality and see at what month does it reach its highest and lowest point. The first one is the city hotel.

city_monthly <- hotel_agg_city %>% 
  filter(market_segment == "Online TA") %>% 
  mutate(date = floor_date(arrival_date, "month")) %>% 
  group_by(date) %>% 
  summarise(demand = sum(demand)) %>% 
  ungroup() 

monthly_ts <- ts(city_monthly,frequency = 12) %>% 
  decompose()

monthly_ts$seasonal %>% 
  matrix(ncol = 12, byrow = T) %>%
  t() %>% 
  as.data.frame() %>% 
  select(V1) %>% 
  mutate(month = month(head(city_monthly$date, 12), label = T)) %>% 
  ggplot(aes(month, V1, fill = V1)) +
  geom_col(show.legend = F) +
  geom_hline(yintercept = 0 ) +
  labs(x = NULL, y = "seasonality", title = "Monthly Seasonality of City Hotel, Online TA Segment ") +
  scale_fill_gradient2(low = "firebrick4", mid = "skyblue", high = "dodgerblue4") +
  theme_pander()

The next one is the online TA segment in Resort Hotel. The seasonality reach its highest point during October and same with the city hotel, it reach the lowest point on March.

resort_monthly <- hotel_agg_resort %>% 
  filter(market_segment == "Online TA") %>% 
  mutate(date = floor_date(arrival_date, "month")) %>% 
  group_by(date) %>% 
  summarise(demand = sum(demand)) %>% 
  ungroup() 

monthly_ts <- ts(resort_monthly,frequency = 12) %>% 
  decompose()

monthly_ts$seasonal %>% 
  matrix(ncol = 12, byrow = T) %>%
  t() %>% 
  as.data.frame() %>% 
  select(V1) %>% 
  mutate(month = month(head(resort_monthly$date, 12), label = T)) %>% 
  ggplot(aes(month, V1, fill = V1)) +
  geom_col(show.legend = F) +
  geom_hline(yintercept = 0 ) +
  labs(x = NULL, y = "seasonality", title = "Monthly Seasonality of resort Hotel, Online TA Segment ") +
  scale_fill_gradient2(low = "firebrick4", mid = "skyblue", high = "dodgerblue4") +
  theme_pander()

Based on both graphs, the high and positive seasonality happens around May-June and September-October. The highest negative seasonality happens in March. Both Lisbon and Algarve are located in Portugal. According to Audley Travels , the best time to visit Portugal is in spring (March-May), when the country is in bloom and waking after the winter. You could also go in fall (between September and October) when the sun is still shining, the weather is warm, and many of the crowds have dispersed. However, the negative seasonality in March and April perhaps tell us that the weather is still too cold to travel around and people love to spend more time to go for vacation during September and October. The summer holiday of school, which is span from late June to early September, have a good influence toward the city hotel seasonality.

September and October are two of the best months to visit Portugal. The weather is still warm and pleasant, and the temperatures are much more manageable for sightseeing or hiking. It’s also a wonderful time to visit many of Portugal’s wineries with the grape harvest in full swing. The beaches are also much quieter.

For the next sections, we will focus on the forecasting of the demand using various machine learning methods.

Preprocessing Specification

We will try several preprocess approach since there is a possibility that transformed data are performing better than the original scale. We will use the following treatment:

• No data transformation
• Squared value
• Scaled value, data will be normalized using z-score
• Log transformation

recipe_spec <- list(
  normal_spec = function(x) x, # no transformation
  squared_spec = function(x) sqrt(x), # square the demand value
  scale_spec = function(x) scale(x), # normalize the demand value with Z-score
  log_spec = function(x) log(x+1) # convert demand to log
) %>% 
  enframe(name = "preprocess_name", value = "preprocess_spec")

recipe_spec

#> # A tibble: 4 × 2
#>   preprocess_name preprocess_spec
#>   <chr>           <list>         
#> 1 normal_spec     <fn>           
#> 2 squared_spec    <fn>           
#> 3 scale_spec      <fn>           
#> 4 log_spec        <fn>

We also need to create a function to reverse the scaled to value into the original value for later model evaluation. For example, if the data is being square rooted, then we need to scale it back by squaring the data.

# reverse function to scale back for model evaluation
reverse_spec <- list(
  
  normal_spec = function(x, y) {
    y <- y
    return(x)
    },
  
  squared_spec = function(x, y) {
    y <- y
    return(x^2)
    },
  
  scale_spec = function(x, y) x * sd(y) + mean(y),
  
  log_spec = function(x,y) {
    y <- y
    return(exp(x)-1)
  }
) %>% 
  enframe(name = "reverse_name", value = "reverse_spec")

# joint the preprocess and the scale-back function
recipe_spec <- recipe_spec %>% 
  bind_cols(reverse_spec)

Seasonality Specification

The next step is to specify the seasonal period for the series. We will try several seasonal period, including:

• weekly seasonality
• monthly seasonality
• annual seasonality
• weekly and monthly seasonality (multi-seasonal)
• weekly and annual seasonality (multi-seasonal)

seasonal_forecast <- list(
  weekly = function(x) ts(x, frequency = 7),
  monthly = function(x) ts(x, frequency = 7*4),
  weekly_monthly = function(x) msts(x, seasonal.periods = c(7, 7*4)),
  weekly_annual = function(x) msts(x, seasonal.periods = c(7, 365)),
  annual = function(x) ts(x, frequency = 365)
  ) %>% 
  enframe(name = "season_name", value = "season_spec")

Impute Outlier

We will also try to preprocess the data by whether an outlier should be replaced or not. If the outlier is replaced, we will identify the outlier and estimate the replacement using the tsoutliers function. Residuals are identified by fitting a loess curve for non-seasonal data and via a periodic STL decomposition for seasonal data.

outlier_spec <- list(
  normal_spec = function(x) x, # no transformation
  
  out_spec = function(x){
    outlier_place <- tsoutliers(x)
    x[outlier_place$index] <- outlier_place$replacement
    return(x)
  }
) %>% 
  enframe(name = "outlier_name", value = "out_spec")

Model Specification

Next, we will specify the model the data will be fit into. The model includes:

• ARIMA
• STL + ETS model
• STL + ARIMA

method_forecast <- list(
  arima  = function(x) auto.arima(x),
  stl_ets = function(x) stlm(x, method = "ets"),
  stl_arima = function(x) stlm(x, method = "arima")
  ) %>% 
  enframe(name = "model_name", value = "model_spec")

Model Fitting

Below is the combination for each specification on each series. For 6 different series, we will have 5 different seasonality specification, 3 different models and other specifications. Therefore, we will 540 different combinations. We will choose the best model based on the RMSE and MAE value on the testing dataset. For example, the first row is the time series for City Hotel with Direct Segment. This data will be transformed using log (log_spec) with annual seasonality and fitted with time series model using ARIMA.

# combine the data with the specification
train_crossing <- crossing(train_list, recipe_spec, seasonal_forecast, outlier_spec, method_forecast )

test_crossing <- crossing(test_list, recipe_spec, seasonal_forecast, outlier_spec, method_forecast)

train_crossing %>% 
  head()

#> # A tibble: 6 × 12
#>   series     data      preprocess_name preprocess_spec reverse_name reverse_spec
#>   <chr>      <list>    <chr>           <list>          <chr>        <list>      
#> 1 City Hote… <tibble … log_spec        <fn>            log_spec     <fn>        
#> 2 City Hote… <tibble … log_spec        <fn>            log_spec     <fn>        
#> 3 City Hote… <tibble … log_spec        <fn>            log_spec     <fn>        
#> 4 City Hote… <tibble … log_spec        <fn>            log_spec     <fn>        
#> 5 City Hote… <tibble … log_spec        <fn>            log_spec     <fn>        
#> 6 City Hote… <tibble … log_spec        <fn>            log_spec     <fn>        
#> # … with 6 more variables: season_name <chr>, season_spec <list>,
#> #   outlier_name <chr>, out_spec <list>, model_name <chr>, model_spec <list>

The following code produce the transformation process for the data before fitted into model.

# model fitting and evaluation
transformed_data <- map2(.x = train_crossing$data,
                         .y = train_crossing$preprocess_spec,
                         .f = ~exec( .y, .x) 
                         )

The following code do all process from transforming data into time series, fitting them into the model and forecast demands for the next 30 days.

# fitting and forecasting
forecast_map <- transformed_data %>% 
  
  # Convert data to time series with different seasonalities
  map2(.y = train_crossing$season_spec,
       .f = ~exec(.y, .x)) %>% 
  
  # Transform outlier
  map2(.y = train_crossing$out_spec,
       .f = ~exec(.y, .x)) %>% 
  
  # Fit data to time series model
  map2(.y = train_crossing$model_spec,
       .f =  ~exec(.y,.x)) %>% 
  
  # Forecast for the next 30 days
  map(forecast, h = 30) %>%  
  
  # Take the mean of the forecast
  map(~pluck(.x, "mean")) %>% 
  
  map(as.numeric)

Forecasting Result

Below is the result for our modeling process for the first time series.

# Result of the first time series
forecast_map[[1]]

#>  [1] 2.346707 2.338587 2.338587 2.338587 2.338587 2.338587 2.338587 2.338587
#>  [9] 2.338587 2.338587 2.338587 2.338587 2.338587 2.338587 2.338587 2.338587
#> [17] 2.338587 2.338587 2.338587 2.338587 2.338587 2.338587 2.338587 2.338587
#> [25] 2.338587 2.338587 2.338587 2.338587 2.338587 2.338587

We use MAE and RMSE to measures and compares the performance of each model.

# scale-back the data to original scale value
forecast_trans <- list()

for (i in 1:length(forecast_map)) {
  forecast_trans[[i]] <- train_crossing$reverse_spec[[i]]( x = forecast_map[[i]], y = transformed_data[[i]])
}

forecast_trans[[1]]

#>  [1] 9.451100 9.366579 9.366579 9.366579 9.366579 9.366579 9.366579 9.366579
#>  [9] 9.366579 9.366579 9.366579 9.366579 9.366579 9.366579 9.366579 9.366579
#> [17] 9.366579 9.366579 9.366579 9.366579 9.366579 9.366579 9.366579 9.366579
#> [25] 9.366579 9.366579 9.366579 9.366579 9.366579 9.366579

# calculate MAE
mae_list <- forecast_trans %>% 
  map2(.y = test_crossing$data, 
       .f = ~yardstick::mae_vec(.x %>% as.numeric(), .y$demand))

rmse_list <- forecast_trans %>% 
  map2(.y = test_crossing$data, 
       .f = ~yardstick::rmse_vec(.x %>% as.numeric(), .y$demand))

# show result
train_crossing %>% 
  separate(series, c("hotel", "market_segment"), sep = "_") %>% 
  select_if(is.character) %>% 
  bind_cols(mae = mae_list %>% as.numeric()) %>% 
  bind_cols(rmse = rmse_list %>% as.numeric()) %>% 
  select(hotel, market_segment, mae, rmse, everything()) %>% 
  head(10)

#> # A tibble: 10 × 9
#>    hotel     market_segment   mae  rmse preprocess_name reverse_name season_name
#>    <chr>     <chr>          <dbl> <dbl> <chr>           <chr>        <chr>      
#>  1 City Hot… Direct          3.79  4.63 log_spec        log_spec     annual     
#>  2 City Hot… Direct         11.7  13.9  log_spec        log_spec     annual     
#>  3 City Hot… Direct         10.6  12.8  log_spec        log_spec     annual     
#>  4 City Hot… Direct          3.79  4.63 log_spec        log_spec     annual     
#>  5 City Hot… Direct         11.7  13.9  log_spec        log_spec     annual     
#>  6 City Hot… Direct         10.6  12.8  log_spec        log_spec     annual     
#>  7 City Hot… Direct          3.79  4.63 log_spec        log_spec     monthly    
#>  8 City Hot… Direct          3.85  4.71 log_spec        log_spec     monthly    
#>  9 City Hot… Direct          3.79  4.64 log_spec        log_spec     monthly    
#> 10 City Hot… Direct          3.79  4.63 log_spec        log_spec     monthly    
#> # … with 2 more variables: outlier_name <chr>, model_name <chr>

Below is the best configuration for each series based on the lowest RMSE, since RMSE give more penality towards large error. To interpret the MAE values, we need to consider the range of the data, shown as the standar deviation of the data.

# best configuration for each series
best_adjust <- train_crossing %>% 
  separate(series, c("hotel", "market_segment"), sep = "_") %>% 
  bind_cols(mae = mae_list %>% as.numeric()) %>%
  bind_cols(rmse = rmse_list %>% as.numeric()) %>% 
  group_by(hotel, market_segment) %>% 
  arrange(rmse) %>% 
  slice(1) 

# show the result
metric_crossing <- train_list %>% crossing(list(mean = mean, std_dev = sd) %>% enframe(name = "type", value = "metric"))

metric_crossing %>% 
  bind_cols(value = map2(.x = metric_crossing$data, 
                         .y = metric_crossing$metric, 
                         .f = ~exec(.y, unlist(.x))) %>% unlist() 
            ) %>% 
  select(series, type, value) %>% 
  pivot_wider(names_from = type, values_from = value) %>% 
  separate(series, c("hotel", "market_segment"), sep = "_") %>% 
  left_join(best_adjust) %>% 
  select_if(~is.list(.) == F) %>% 
  select(-reverse_name) %>% 
  select(hotel, market_segment, mae, rmse,  mean, std_dev, everything()) %>% 
  arrange(rmse)

#> # A tibble: 6 × 10
#>   hotel     market_segment   mae  rmse  mean std_dev preprocess_name season_name
#>   <chr>     <chr>          <dbl> <dbl> <dbl>   <dbl> <chr>           <chr>      
#> 1 Resort H… Direct          2.21  2.65  8.14    4.85 normal_spec     weekly     
#> 2 City Hot… Direct          3.43  4.17  7.53    4.83 normal_spec     weekly     
#> 3 Resort H… Offline TA/TO   4.26  6.09  9.39    7.38 log_spec        weekly_ann…
#> 4 Resort H… Online TA       8.71 11.0  21.8    11.7  squared_spec    weekly     
#> 5 City Hot… Offline TA/TO   8.80 13.7  21.4    32.0  normal_spec     weekly     
#> 6 City Hot… Online TA      13.2  16.2  48.1    28.5  squared_spec    weekly     
#> # … with 2 more variables: outlier_name <chr>, model_name <chr>

Since we don’t have much data, only two years of transactions, the model performance may not perform so well. However, judging from the MAE value, the performance is quite acceptable, with most of the error values are less than the value of one standard deviation. Compared to the aggregated data on the first forecast which have MAE of 21, we have lower MAE value for each series, which give us an evidence that by making a separate forecasting models for each market segment will make the model more accurate.

# Aggregated time series performance
accuracy(forecast_agg, test_agg$demand) %>%
  as.data.frame() %>% 
  slice(2) %>% 
  select(MAE, RMSE)

#>               MAE     RMSE
#> Test set 21.25028 26.31801

Below is the forecasting result for each series. The red line indicate the actual demand value while the blue line indicate the forecast value. The blue area represent area with 80% prediction interval while the light blue are represent the 95% prediction interval. Most of the actual demand is still inside the forecasting intervals.

model_best <- map2(.x = best_adjust$data, 
                   .y = best_adjust$season_spec,
                   .f =  ~exec(.y,.x)) %>% 
  map2(.y = best_adjust$model_spec,
       .f =  ~exec(.y,.x)
       )

data_test <- test_list$data[[6]] %>% 
  msts(start = 110, seasonal.periods =  7)

model_best[[6]] %>%
  forecast(h = 30 ) %>%  
  autoplot() +
  autolayer(data_test, series = "Data Test") +
  scale_color_manual(values = "firebrick") +
  labs(subtitle = "Resort Hotel, Online TA",
       y = "Demand", x  = NULL) +
  theme_pander() +
  scale_x_continuous(limits = c(100, 115) ,
                     labels = as.Date.numeric(seq(100,115,5)*7-7, origin = range(hotel$arrival_date)[1])) +
  theme(legend.position = "top")

model_best[[6]]$residuals %>% 
  as.data.frame() %>% 
  ggplot(aes(x)) +
  geom_density(fill = "skyblue", alpha = 0.7, color = "white") +
  labs(x = "Residuals", y = "Density",
       title = "Residuals Distribution") +
  theme_pander()

data_test <- test_list$data[[3]] %>% 
  ts(start = 110, frequency = 7)

model_best[[3]] %>%
  forecast(h = 30 ) %>% 
  autoplot() +
  autolayer(data_test, series = "Data Test") +
  scale_color_manual(values = "firebrick") +
  labs(subtitle = "City Hotel, Online TA",
       y = "Demand", x  = NULL) +
  theme_pander() +
  scale_x_continuous(limits = c(100, 115), 
                     labels = as.Date.numeric(seq(100,115,5)*7-7, origin = range(hotel$arrival_date)[1])) +
  theme(legend.position = "top")

model_best[[3]]$residuals %>% 
  as.data.frame() %>% 
  ggplot(aes(x)) +
  geom_density(fill = "skyblue", alpha = 0.7, color = "white") +
  labs(x = "Residuals", y = "Density",
       title = "Residuals Distribution") +
  theme_pander()

data_test <- test_list$data[[5]] %>% msts(start = 3.090411, seasonal.periods = c(7, 365))

model_best[[5]] %>%
  forecast(h = 30 ) %>% 
  autoplot() +
  autolayer(data_test, series = "Data Test") +
  scale_color_manual(values = "firebrick") +
  labs(subtitle = "Resort Hotel, Offline TA/TO",
       y = "Demand", x  = NULL) +
  theme_pander() +
  scale_x_continuous(limits = c(3, 3.2), 
                     labels = as.Date.numeric(seq(3,3.2,0.05)*365-365, origin = range(hotel$arrival_date)[1])) +
  theme(legend.position = "top")

model_best[[5]]$residuals %>% 
  as.data.frame() %>% 
  ggplot(aes(x)) +
  geom_density(fill = "skyblue", alpha = 0.7, color = "white") +
  labs(x = "Residuals", y = "Density",
       title = "Residuals Distribution") +
  theme_pander()

data_test <- test_list$data[[2]] %>% ts(start = 110, frequency = 7)

model_best[[2]] %>%
  forecast(h = 30 ) %>% 
  autoplot() +
  autolayer(data_test, series = "Data Test") +
  scale_color_manual(values = "firebrick") +
  labs(subtitle = "City Hotel, Offline TA/TO",
       y = "Demand", x  = NULL) +
  theme_pander() +
  scale_x_continuous(limits = c(100, 115), 
                     labels = as.Date.numeric(seq(100,115,5)*7-7, origin = range(hotel$arrival_date)[1])) +
  theme(legend.position = "top")

model_best[[2]]$residuals %>% 
  as.data.frame() %>% 
  ggplot(aes(x)) +
  geom_density(fill = "skyblue", alpha = 0.7, color = "white") +
  labs(x = "Residuals", y = "Density",
       title = "Residuals Distribution") +
  theme_pander()

data_test <- test_list$data[[4]] %>% ts(start = 110, frequency = 7)

model_best[[4]] %>%
  forecast(h = 30 ) %>% 
  autoplot() +
  autolayer(data_test, series = "Data Test") +
  scale_color_manual(values = "firebrick") +
  labs(subtitle = "Resort Hotel, Direct",
       y = "Demand", x  = NULL) +
  theme_pander() +
  scale_x_continuous(limits = c(100, 115), 
                     labels = as.Date.numeric(seq(100,115,5)*7-7, origin = range(hotel$arrival_date)[1])) +
  theme(legend.position = "top")

model_best[[4]]$residuals %>% 
  as.data.frame() %>% 
  ggplot(aes(x)) +
  geom_density(fill = "skyblue", alpha = 0.7, color = "white") +
  labs(x = "Residuals", y = "Density",
       title = "Residuals Distribution") +
  theme_pander()

data_test <- test_list$data[[1]] %>% ts(start = 110, frequency = 7)

model_best[[1]] %>%
  forecast(h = 30 ) %>% 
  autoplot() +
  autolayer(data_test, series = "Data Test") +
  scale_color_manual(values = "firebrick") +
  labs(subtitle = "City Hotel, Direct",
       y = "Demand", x  = NULL) +
  theme_pander() +
  scale_x_continuous(limits = c(100, 115), 
                     labels = as.Date.numeric(seq(100,115,5)*7-7, origin = range(hotel$arrival_date)[1])) +
  theme(legend.position = "top")

model_best[[1]]$residuals %>% 
  as.data.frame() %>% 
  ggplot(aes(x)) +
  geom_density(fill = "skyblue", alpha = 0.7, color = "white") +
  labs(x = "Residuals", y = "Density",
       title = "Residuals Distribution") +
  theme_pander()

Autocorrelation

The autocorrelation can be checked using the Ljung-Box test. If there are correlations between residuals, then there is information left in the residuals which should be used in computing forecasts.

best_adjust %>% 
  select(hotel, market_segment) %>% 
  bind_cols( 
      map(model_best, ~Box.test(.x$residuals, type = "Lj")) %>% 
  unlist() %>% 
  matrix(ncol = 5, byrow = T) %>% 
  as.data.frame() %>% 
  select(4:3) %>% 
  rename(p_value = V3, test = V4) %>% 
  mutate(p_value = p_value %>% 
           as.character() %>% 
           as.numeric() %>% 
           round(4)) 
  )

#> # A tibble: 6 × 4
#> # Groups:   hotel, market_segment [6]
#>   hotel        market_segment test           p_value
#>   <chr>        <chr>          <chr>            <dbl>
#> 1 City Hotel   Direct         Box-Ljung test   0.976
#> 2 City Hotel   Offline TA/TO  Box-Ljung test   0.971
#> 3 City Hotel   Online TA      Box-Ljung test   0.769
#> 4 Resort Hotel Direct         Box-Ljung test   0.854
#> 5 Resort Hotel Offline TA/TO  Box-Ljung test   0.959
#> 6 Resort Hotel Online TA      Box-Ljung test   0.959

The results suggests that all of our models don’t have any autocorrelation based on the non-significant p-value.

Normality

We will also check if the residuals for each model is normally distributed using Shapiro-Wilk Test. If the residuals are not normally distributed, it will lead to a biased parameter and less optimal forecast. This is also indicate that we can still tweak our model to get a better performance.

best_adjust %>% 
  select(hotel, market_segment) %>% 
  bind_cols(
    mean_error = map(model_best, ~mean(.$residuals)) %>% 
      as.numeric() %>% 
      round(5)
  ) %>%
  bind_cols(
    median_error = map(model_best, ~median(.$residuals)) %>% 
      as.numeric() %>% 
      round(5)
  ) %>% 
  bind_cols( 
      map(model_best, ~shapiro.test(.x$residuals)) %>% 
        unlist() %>% 
        matrix(ncol = 4, byrow = T) %>% 
        as.data.frame() %>% 
        select(3:2) %>% 
        rename(p_value = V2, test = V3)
  ) %>% 
  rename(segment = market_segment) %>% 
  mutate(p_value = p_value %>% 
           as.character() %>% 
           as.numeric() %>% 
           number(accuracy = 0.00001))

#> # A tibble: 6 × 6
#> # Groups:   hotel, segment [6]
#>   hotel       segment      mean_error median_error test                  p_value
#>   <chr>       <chr>             <dbl>        <dbl> <chr>                 <chr>  
#> 1 City Hotel  Direct           0.180        -0.387 Shapiro-Wilk normali… 0.00000
#> 2 City Hotel  Offline TA/…     0.104        -6.74  Shapiro-Wilk normali… 0.00000
#> 3 City Hotel  Online TA        0.0146       -0.898 Shapiro-Wilk normali… 0.00000
#> 4 Resort Hot… Direct           0.0442       -0.586 Shapiro-Wilk normali… 0.00000
#> 5 Resort Hot… Offline TA/…    -0.0813       -0.142 Shapiro-Wilk normali… 0.00000
#> 6 Resort Hot… Online TA        0.211        -0.426 Shapiro-Wilk normali… 0.00000

Based on the result, all of our model didn’t fulfill the normality assumption for the residuals. The positive mean of error signify that the model is underestimate the forecast while negative mean error means the model is overestimate. If we look at the median of error, all of our models are underestimate on the training set. They might be influenced by the presence of an outlier value such as a really high demand especially on the early part of the series. This suggest that we can improve the model further in order to get better performance.