<jupyter_start><jupyter_text>Data Prep  Production<jupyter_code>production = read_excel('../data/rice_production.xlsx', 'all')

production_annual = production %>% 
  group_by(year, rice_type) %>%
  summarize(area_plant = sum(area_plant),
            area_harvest = sum(area_harvest),
            production = sum(production)*0.66 # convert to milled rice
           ) %>%
  mutate(yield_plant = production*1000/area_plant,
         yield_harvest = production*1000/area_harvest,
        ) %>% ungroup()

production_wide = production_annual %>%
  pivot_wider(
    year,
    rice_type,
    values_from = production
  )

names(production_wide) = c(
  'year_th', 'q_pathum', 'q_hommali', 'q_white', 'q_glutinous'
)

production_wide %>% tail()<jupyter_output>`summarise()` regrouping output by 'year' (override with `.groups` argument)<jupyter_text>Export<jupyter_code>export = read_rds('../data/export_1006_master.rds')
export_meta = read_excel('../data/meta_hs1006_th.xlsx')
export_joined = export %>% left_join(export_meta, by = c('hscode' = 'hscode'))

export_annual = export_joined %>% 
  group_by(year_th, varities) %>%
  summarize(vol = sum(vol)) %>% ungroup()

export_wide = export_annual %>%
  pivot_wider(
    year_th,
    varities,
    values_from = vol
  )

names(export_wide) = c(
  'year_th', 'ex_glutinous', 'ex_white', 
  'ex_hommali', 'ex_pathum', 'ex_colored'
)

head(export_wide)

export_monthly = export_joined %>% 
  group_by(year_th, month, varities) %>%
  summarize(vol = sum(vol)) %>% ungroup()

export_wide_m = export_monthly %>%
  pivot_wider(
    c(year_th, month),
    varities,
    values_from = vol,
    values_fill = 0
  ) %>% filter(year_th >= 2549)

names(export_wide_m) = c(
  'year_th', 'month',
  'ex_glutinous','ex_white','ex_hommali','ex_pathum','ex_colored'
)<jupyter_output>`summarise()` regrouping output by 'year_th' (override with `.groups` argument)<jupyter_text>Domestic Consumption<jupyter_code>data_m = export_wide_m %>% 
  mutate(date = as.Date(paste0(year_th-543, '-', month, '-1')))

data_a = export_wide %>% 
  left_join(production_wide, by = 'year_th') %>%
  mutate(dom_hommali = q_hommali - ex_hommali,
         dom_pathum = q_pathum - ex_pathum,
         dom_white = q_white - ex_white,
         dom_glutinous = q_glutinous - ex_glutinous
        ) %>%
  filter(year_th >= 2549)

# save in csv format
data_a %>% write_excel_csv('../data/rice_a.csv')
data_m %>% write_excel_csv('../data/rice_m.csv')
# export_names = c('ex_glutinous', 'ex_white', 'ex_hommali', 'ex_pathum')
# ts_export = ts(data_m[export_names], start = 2549-543, frequency = 12)<jupyter_output><empty_output><jupyter_text>Read<jupyter_code>rice_a = read_csv('../data/rice_a.csv') 
rice_m = read_csv('../data/rice_m.csv')
rice_m_ts = rice_m %>% tk_ts(start = 2549-543, frequency = 12)<jupyter_output>Parsed with column specification:
cols(
  year_th = [32mcol_double()[39m,
  ex_glutinous = [32mcol_double()[39m,
  ex_white = [32mcol_double()[39m,
  ex_hommali = [32mcol_double()[39m,
  ex_pathum = [32mcol_double()[39m,
  ex_colored = [32mcol_double()[39m,
  q_pathum = [32mcol_double()[39m,
  q_hommali = [32mcol_double()[39m,
  q_white = [32mcol_double()[39m,
  q_glutinous = [32mcol_double()[39m,
  dom_hommali = [32mcol_double()[39m,
  dom_pathum = [32mcol_double()[39m,
  dom_white = [32mcol_double()[39m,
  dom_glutinous = [32mcol_double()[39m
)

Parsed with column specification:
cols(
  year_th = [32mcol_double()[39m,
  month = [31mcol_character()[39m,
  ex_glutinous = [32mcol_double()[39m,
  ex_white = [32mcol_double()[39m,
  ex_hommali = [32mcol_double()[39m,
  ex_pathum = [32mcol_double()[39m,
  ex_colored = [32mcol_double()[39m,
  date = [34mcol_date(format = "")[39m
)

Warning message:
"Non-numeric columns being dropped: month, date"<jupyter_text>Data Visualization<jupyter_code>p_hommali = ggplot(rice_a, aes(x=year_th)) +
  geom_line( aes(y=dom_hommali/1e6, color = 'Domestic')) + 
  geom_line( aes(y=ex_hommali/1e6, color = 'Export')) + 
  scale_colour_manual("",
                      label = c('บริโภคในประเทศ', 'ส่งออก'),
                      breaks = c("Domestic", "Export"),
                      values = c("Domestic"="#0b53c1", 
                                 "Export"="#ff0055")
                     ) +
  labs(title = 'ข้าวหอมมะลิ',
       x =  "ปี",
       y = 'ล้านตันข้าวสาร'
      ) +
  expand_limits(y = 0) +
  scale_y_continuous(label = comma) +
  theme_ipsum(
    base_size = 16,
    base_family = "Athiti Light",
    axis_title_size = 16
  ) +
  theme(
    legend.position = c(1, 1),
    legend.direction = 'horizontal',
    legend.justification = 'right',
    legend.title = element_blank()
  )

p_hommali
ggsave('../figures/plot_hommali.png', width = 8, height = 4)
p_pathum = ggplot(rice_a, aes(x=year_th)) +
  geom_line( aes(y=dom_pathum/1e6, color = 'Domestic')) + 
  geom_line( aes(y=ex_pathum/1e6, color = 'Export')) + 
  scale_colour_manual("",
                      label = c('บริโภคในประเทศ', 'ส่งออก'),
                      breaks = c("Domestic", "Export"),
                      values = c("Domestic"="#0b53c1", 
                                 "Export"="#ff0055")
                     ) +
  labs(title = 'ข้าวหอมปทุม',
       x =  "ปี",
       y = 'ล้านตันข้าวสาร'
      ) +
  expand_limits(y = 0) +
  scale_y_continuous(label = comma) +
  theme_ipsum(
    base_size = 16,
    base_family = "Athiti Light",
    axis_title_size = 16
  ) +
  theme(
    legend.position = c(1, 1),
    legend.direction = 'horizontal',
    legend.justification = 'right',
    legend.title = element_blank()
  )

p_pathum
ggsave('../figures/plot_pathum.png', width = 8, height = 4)
p_white = ggplot(rice_a, aes(x=year_th)) +
  geom_line( aes(y=dom_white/1e6, color = 'Domestic')) + 
  geom_line( aes(y=ex_white/1e6, color = 'Export')) + 
  scale_colour_manual("",
                      label = c('บริโภคในประเทศ', 'ส่งออก'),
                      breaks = c("Domestic", "Export"),
                      values = c("Domestic"="#0b53c1", 
                                 "Export"="#ff0055")
                     ) +
  labs(title = 'ข้าวเจ้าขาว',
       x =  "ปี",
       y = 'ล้านตันข้าวสาร'
      ) +
  expand_limits(y = 0) +
  scale_y_continuous(label = comma) +
  theme_ipsum(
    base_size = 16,
    base_family = "Athiti Light",
    axis_title_size = 16
  ) +
  theme(
    legend.position = c(1, 1),
    legend.direction = 'horizontal',
    legend.justification = 'right',
    legend.title = element_blank()
  )

p_white
ggsave('../figures/plot_white.png', width = 8, height = 4)
p_glutinous = ggplot(rice_a, aes(x=year_th)) +
  geom_line( aes(y=dom_glutinous/1e6, color = 'Domestic')) + 
  geom_line( aes(y=ex_glutinous/1e6, color = 'Export')) + 
  scale_colour_manual("",
                      label = c('บริโภคในประเทศ', 'ส่งออก'),
                      breaks = c("Domestic", "Export"),
                      values = c("Domestic"="#0b53c1", 
                                 "Export"="#ff0055")
                     ) +
  labs(title = 'ข้าวเหนียว',
       x =  "ปี",
       y = 'ล้านตันข้าวสาร'
      ) +
  expand_limits(y = 0) +
  scale_y_continuous(label = comma) +
  theme_ipsum(
    base_size = 16,
    base_family = "Athiti Light",
    axis_title_size = 16
  ) +
  theme(
    legend.position = c(1, 1),
    legend.direction = 'horizontal',
    legend.justification = 'right',
    legend.title = element_blank()
  )

p_glutinous
ggsave('../figures/plot_glutinous.png', width = 8, height = 4)
ggplot(rice_m, aes(x=date)) +
  geom_line( aes(y=ex_hommali/1e6, color = 'ข้าวหอมมะลิ')) +
  geom_line( aes(y=ex_pathum/1e6, color = 'ข้าวหอมปทุม')) +
  geom_line( aes(y=ex_white/1e6, color = 'ข้าวเจ้าขาว')) +
  geom_line( aes(y=ex_glutinous/1e6, color = 'ข้าวเหนียว')) +
  scale_color_d3() +
  expand_limits(y = 0) +
  scale_y_continuous(label = comma) +
  labs(title = 'ปริมาณส่งออกข้าวสารของไทย',
       x =  "ปี",
       y = 'ล้านตันข้าวสาร'
      ) +
  theme_ipsum(
    base_size = 16,
    base_family = "Athiti Light",
    axis_title_size = 16
  ) +
  theme(
    legend.position = c(1, 1),
    legend.direction = 'horizontal',
    legend.justification = 'right',
    legend.title = element_blank()
  ) +
  scale_x_date(date_breaks = '2 years',
               labels = function(x) year(x)+543
               )

ggsave('../figures/plot_export_m.png', width = 14, height = 8)<jupyter_output><empty_output><jupyter_text>Unit Root Test<jupyter_code>dom_names = c('dom_hommali', 'dom_pathum', 'dom_white', 'dom_glutinous')
data_dom = ts(data[, dom_names])
get_objname = function(x) deparse(substitute(x))

ur_adf = function(y, n_diff = 0, varname = NULL, ...) {
  
  varname = if (is.null(varname)) deparse(substitute(y)) else varname
  y = if(n_diff == 0) y else diff(y, n_diff)
  
  ur.trend = ur.df(y, type='trend', selectlags = "AIC", ...)
  ur.drift = ur.df(y, type='drift', selectlags = "AIC", ...)
  ur.none  = ur.df(y, type='none', selectlags = "AIC", ...)

  tstat.trend = ur.trend@teststat
  tstat.drift = ur.drift@teststat
  tstat.none  = ur.none@teststat

  cv.trend = ur.trend@cval
  cv.drift = ur.drift@cval
  cv.none  = ur.none@cval

  df_test = rbind(
    cbind(t(tstat.trend), cv.trend),
    cbind(t(tstat.drift), cv.drift),
    cbind(t(tstat.none) , cv.none)
  ) %>% 
  as.data.frame() %>%
  rownames_to_column("hypo") %>%
  mutate(
    result = ifelse(abs(statistic) >= abs(`5pct`), 'Reject', 'Accept'),
    variable = varname,
    level = paste0('d', n_diff),
    statistic = round(statistic, 2),
  ) %>%
  filter(str_starts(hypo, 'tau')) %>%
  select(variable, hypo, level, everything())
  
  return(df_test)
}
dom_hommali   = rice_a$dom_hommali
dom_pathum    = rice_a$dom_pathum
dom_white     = rice_a$dom_white
dom_glutinous = rice_a$dom_glutinous

ex_hommali    = rice_m_ts[, 'ex_hommali']
ex_pathum     = rice_m_ts[, 'ex_pathum'] 
ex_white      = rice_m_ts[, 'ex_white']
ex_glutinous  = rice_m_ts[, 'ex_glutinous']

rbind(ur_adf(dom_hommali, 0)  , ur_adf(dom_hommali, 1),
      ur_adf(dom_pathum, 0)   , ur_adf(dom_pathum, 1),
      ur_adf(dom_white, 0)    , ur_adf(dom_white, 1),
      ur_adf(dom_glutinous, 0), ur_adf(dom_glutinous, 1),
      ur_adf(ex_hommali, 0)   , ur_adf(ex_hommali, 1),
      ur_adf(ex_pathum, 0)    , ur_adf(ex_pathum, 1),
      ur_adf(ex_white, 0)     , ur_adf(ex_white, 1),
      ur_adf(ex_glutinous, 0) , ur_adf(ex_glutinous, 1)
     ) %>% print()

# write_excel_csv('../results/adf_table.csv')<jupyter_output>variable hypo level statistic  1pct  5pct 10pct result
1    dom_hommali tau3    d0     -1.83 -4.38 -3.60 -3.24 Accept
2    dom_hommali tau2    d0     -1.74 -3.75 -3.00 -2.63 Accept
3    dom_hommali tau1    d0      1.19 -2.66 -1.95 -1.60 Accept
4    dom_hommali tau3    d1     -1.75 -4.38 -3.60 -3.24 Accept
5    dom_hommali tau2    d1     -1.78 -3.75 -3.00 -2.63 Accept
6    dom_hommali tau1    d1     -1.41 -2.66 -1.95 -1.60 Accept
7     dom_pathum tau3    d0     -4.27 -4.38 -3.60 -3.24 Reject
8     dom_pathum tau2    d0     -1.36 -3.75 -3.00 -2.63 Accept
9     dom_pathum tau1    d0     -0.66 -2.66 -1.95 -1.60 Accept
10    dom_pathum tau3    d1     -1.83 -4.38 -3.60 -3.24 Accept
11    dom_pathum tau2    d1     -3.63 -3.75 -3.00 -2.63 Reject
12    dom_pathum tau1    d1     -3.10 -2.66 -1.95 -1.60 Reject
13     dom_white tau3    d0     -2.60 -4.38 -3.60 -3.24 Accept
14     dom_white tau2    d0     -2.34 -3.75 -3.00 -2.63 Accept
15     dom_white tau1    d0     -0.87 -2.66 -1.95 -1.60[...]<jupyter_text>Define function<jupyter_code># fit function
fets = function(x, model='AAN', damped = FALSE, h=h) {
  forecast::ets(x, model = model, damped = damped) %>% forecast(h=h)
}

fets.auto = function(x, h=h) {
  forecast::ets(x, model = 'ZZZ', damped = NULL) %>% forecast(h=h)
}

farima = function(x, order = c(0,0,0), seasonal = c(0,0,0), h=h) {
  forecast::Arima(x,
                  order = order, 
                  seasonal = seasonal,
                 ) %>% forecast(h=h)
}

farima.auto = function(x, d=0, h=h) {
  forecast::auto.arima(x, d=d) %>% forecast(h=h)
}


# accuracy measurement
rmse = function(data, h = h) {
  data = data %>% select('actual', starts_with('yhat')) %>% exp()
  h = ncol(data) - 1
  for (i in seq(h)) {
    data[, paste0('e', i)] = (
      lead(data[,'actual'], i) - data[, paste0('yhat', i)]
    )^2
  }
  rmse_score = data %>% select(starts_with('e')) %>% 
    as.matrix() %>% mean(na.rm = TRUE) %>% sqrt()
  
  return(rmse_score)
}

mape = function(data, h = h) {
  data = data %>% select('actual', starts_with('yhat')) %>% exp()
  h = ncol(data) - 1
  for (i in seq(h)) {
    data[, paste0('e', i)] = abs(
      lead(data[,'actual'], i) - data[, paste0('yhat', i)]
    )*100 / lead(data[,'actual'], i)
  }
  mape_score = data %>% select(starts_with('e')) %>% 
    as.matrix() %>% mean(na.rm = TRUE)
  
  return(mape_score)
}

eval_ets = function(x, params, h = 2, initial = 11) {
  pmap(
    list(params$model, params$damped), 
    function(first, second) {
      
      # 1 compute forecast error using time series cross validation
      error = tsCV(x, fets, h = h, initial = initial, 
                   model = first, 
                   damped = second
                  ) %>% as.data.frame()
      new_colnames = c('actual', paste0('e', seq(1:h)))
      error = data.frame(x, error) %>% set_names(new_colnames)
      
      # 2 compute predicted (yhat)
      for (i in seq(h)) {
        error[, paste0('yhat', i)] = (
          lead(error[,'actual'], i) - error[, paste0('e', i)]
        )
      }
      
      # 3 compute rmse, mape
      data.frame(
        model = paste(first, second, sep = '-'),
        rmse = rmse(error),
        mape = mape(error)
      ) %>% na.omit()
    }
  ) %>% 
  reduce(bind_rows) %>%
  arrange(rmse)
}

eval_arima = function(x, params, h = 2, initial = 11) {
  pmap(
    list(params$order, params$seasonal), 
    function(first, second) {
      
      # 1 compute forecast error using time series cross validation
      error = tsCV(x, farima, h = h, initial = initial, 
                   order = first, 
                   seasonal = second
                  ) %>% as.data.frame()
      new_colnames = c('actual', paste0('e', seq(1:h)))
      error = data.frame(x, error) %>% set_names(new_colnames)
      
      # 2 compute predicted (yhat)
      for (i in seq(h)) {
        error[, paste0('yhat', i)] = (
          lead(error[,'actual'], i) - error[, paste0('e', i)]
        )
      }
  
      # 3 compute rmse, mape
      data.frame(
        model = paste(paste(first, collapse = ','), 
                      paste(second, collapse = ','), sep = '-'),
        rmse = rmse(error),
        mape = mape(error)
      ) %>% na.omit()
    }
  ) %>% 
  reduce(bind_rows) %>%
  arrange(rmse)
}<jupyter_output><empty_output><jupyter_text>Domestic<jupyter_code>Z1 = c('A','M')
Z2 = c('N', 'A', 'M')
Z3 = c('N')
damped = c(FALSE, TRUE)

ets_params = expand.grid(Z1, Z2, Z3, damped) %>%
  mutate(
    model = paste0(Var1, Var2, Var3),
    damped = Var4
) %>% select(model, damped)
length(dom_hommali) * 0.8
e.ets.dom_hml = eval_ets(
  dom_hommali %>% log(), params = ets_params, h = 2, initial = 11
)
e.ets.dom_pat = eval_ets(
  dom_pathum %>% log(), params = ets_params, h = 2, initial = 11
)
e.ets.dom_wht = eval_ets(
  dom_white %>% log(), params = ets_params, h = 2, initial = 11
)
e.ets.dom_glu = eval_ets(
  dom_glutinous %>% log(), params = ets_params, h = 2, initial = 11
)<jupyter_output><empty_output><jupyter_text>Export<jupyter_code># ETS params
Z1 = c('A','M')
Z2 = c('N', 'A', 'M')
Z3 = c('N', 'A', 'M')
damped = c(FALSE, TRUE)

ets_params = expand.grid(Z1, Z2, Z3, damped) %>%
  mutate(
    model = paste0(Var1, Var2, Var3),
    damped = Var4
) %>% select(model, damped)


# ARIMA params
p = seq(0,2)
q = seq(0,2)
P = seq(0,1)
Q = seq(0,1)

arima_params = expand.grid(p,q,P,Q) %>% rowwise() %>%
  mutate(
    order = pmap(list(Var1, Var2), function(x, y) c(x, 0, y)),
    seasonal = pmap(list(Var3, Var4), function(x, y) c(x, 0, y)),
  ) %>% select(order, seasonal)
length(ex_hommali) * 0.8 / 12
11*12
h = 24
initial = 132

e.ets.ex_hml = eval_ets(
  ex_hommali %>% log(), params = ets_params, h = h, initial = initial
)
e.ets.ex_pat = eval_ets(
  ex_pathum %>% log(), params = ets_params, h = h, initial = initial
)
e.ets.ex_wht = eval_ets(
  ex_white %>% log(), params = ets_params, h = h, initial = initial
)
e.ets.ex_glu = eval_ets(
  ex_glutinous %>% log(), params = ets_params, h = h, initial = initial
)
e.arima.ex_hml = eval_arima(
  ex_hommali %>% log(), params = arima_params, h = h, initial = initial
)
e.arima.ex_pat = eval_arima(
  ex_pathum %>% log(), params = arima_params, h = h, initial = initial
)
e.arima.ex_wht = eval_arima(
  ex_white %>% log(), params = arima_params, h = h, initial = initial
)
e.arima.ex_glu = eval_arima(
  ex_glutinous %>% log(), params = arima_params, h = h, initial = initial
)<jupyter_output><empty_output><jupyter_text>Domestic Consumption<jupyter_code>x = dom_hommali
e.ets = tsCV(log(x), fets, h = 2, initial = 10)
ets.RMSE = sqrt(mean(e.ets^2, na.rm=TRUE))
# e.arima = tsCV(log(x), farima, d=1, h = 2, initial = 10)
# arima.RMSE = sqrt(mean(e.arima^2, na.rm=TRUE))

ets.RMSE
# arima.RMSE
x = dom_hommali
e.ets = tsCV(log(x), fets, h = 2, initial = 10)
ets.RMSE = sqrt(mean(e.ets^2, na.rm=TRUE))
e.arima = tsCV(log(x), farima, d=1, h = 2, initial = 10)
arima.RMSE = sqrt(mean(e.arima^2, na.rm=TRUE))

ets.RMSE
arima.RMSE
x = data$dom_pathum
e.ets = tsCV(log(x), fets, h = 2, initial = 10)
ets.RMSE = sqrt(mean(e.ets^2, na.rm=TRUE))
e.arima = tsCV(log(x), farima, d=1, h = 2, initial = 10)
arima.RMSE = sqrt(mean(e.arima^2, na.rm=TRUE))

ets.RMSE
arima.RMSE
x = data$dom_white
e.ets = tsCV(log(x), fets, h = 2, initial = 10)
ets.RMSE = sqrt(mean(e.ets^2, na.rm=TRUE))
e.arima = tsCV(log(x), farima, d=1, h = 2, initial = 10)
arima.RMSE = sqrt(mean(e.arima^2, na.rm=TRUE))

ets.RMSE
arima.RMSE
x = data$dom_glutinous
e.ets = tsCV(log(x), fets, h = 2, initial = 10)
ets.RMSE = sqrt(mean(e.ets^2, na.rm=TRUE))
e.arima = tsCV(log(x), farima, d=1, h = 2, initial = 10)
arima.RMSE = sqrt(mean(e.arima^2, na.rm=TRUE))

ets.RMSE
arima.RMSE<jupyter_output><empty_output><jupyter_text>Export<jupyter_code>x = ts_export[, 'ex_hommali']
e.ets = tsCV(log(x), fets, h = 24, initial = 120)
ets.RMSE = sqrt(mean(e.ets^2, na.rm=TRUE))
e.arima = tsCV(log(x), farima, d=1, h = 24, initial = 120)
arima.RMSE = sqrt(mean(e.arima^2, na.rm=TRUE))

ets.RMSE
arima.RMSE
x = ts_export[, 'ex_pathum']
e.ets = tsCV(log(x), fets, h = 24, initial = 120)
ets.RMSE = sqrt(mean(e.ets^2, na.rm=TRUE))
e.arima = tsCV(log(x), farima, d=1, h = 24, initial = 120)
arima.RMSE = sqrt(mean(e.arima^2, na.rm=TRUE))

ets.RMSE
arima.RMSE
x = ts_export[, 'ex_white']
e.ets = tsCV(log(x), fets, h = 24, initial = 120)
ets.RMSE = sqrt(mean(e.ets^2, na.rm=TRUE))
e.arima = tsCV(log(x), farima, d=1, h = 24, initial = 120)
arima.RMSE = sqrt(mean(e.arima^2, na.rm=TRUE))

ets.RMSE
arima.RMSE
x = ts_export[, 'ex_glutinous']
e.ets = tsCV(log(x), fets, h = 24, initial = 120)
ets.RMSE = sqrt(mean(e.ets^2, na.rm=TRUE))
e.arima = tsCV(log(x), farima, d=1, h = 24, initial = 120)
arima.RMSE = sqrt(mean(e.arima^2, na.rm=TRUE))

ets.RMSE
arima.RMSE
arima
fets <- function(x, model='AAN', damped = FALSE, h=1) {
  forecast(ets(x, model=model, damped = damped), h=h)
}

farima <- function(x, d, h=1) {
  forecast(auto.arima(x, d=d), h=h)
}
fets <- function(x, model='ZZZ', damped = NULL, h=2) {
  forecast(ets(x, model=model, damped = damped), h=h)
}
fets(log(data$dom_hommali)) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Hommali Rice')
ggsave('../figures/plot_ets_hommali_d.png', width = 8, height = 4)

fets(log(data$dom_pathum)) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Pathumtani Rice')
ggsave('../figures/plot_ets_pathum_d.png', width = 8, height = 4)

fets(log(data$dom_white)) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'White Rice')
ggsave('../figures/plot_ets_white_d.png', width = 8, height = 4)

fets(log(data$dom_glutinous)) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Glutinous Rice')
ggsave('../figures/plot_ets_glutinous_d.png', width = 8, height = 4)
fets(log(ts_export[,'ex_hommali']), h = 24) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Hommali Rice')
ggsave('../figures/plot_ets_exp_hml.png', width = 8, height = 4)

fets(log(ts_export[,'ex_pathum']), h = 24) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Pathumtani Rice')
ggsave('../figures/plot_ets_exp_pat.png', width = 8, height = 4)

fets(log(ts_export[,'ex_white']), h = 24) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'White Rice')
ggsave('../figures/plot_ets_exp_wht.png', width = 8, height = 4)

fets(log(ts_export[,'ex_glutinous']), h = 24) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Glutinous Rice')
ggsave('../figures/plot_ets_exp_glu.png', width = 8, height = 4)<jupyter_output><empty_output><jupyter_text>ARIMA<jupyter_code>farima <- function(x, d, h=2) {
  forecast(auto.arima(x, d=d), h=h)
}
farima(log(data$dom_hommali), d=1, h=2) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Hommali Rice')
ggsave('../figures/plot_arima_dom_hml.png', width = 8, height = 4)

farima(log(data$dom_pathum), d=1, h=2) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Pathumtani Rice')
ggsave('../figures/plot_arima_dom_pat.png', width = 8, height = 4)

farima(log(data$dom_white), d=1, h=2) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'White Rice')
ggsave('../figures/plot_arima_dom_wht.png', width = 8, height = 4)

farima(log(data$dom_glutinous), d=1, h=2) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Glutinous Rice')
ggsave('../figures/plot_arima_dom_glu.png', width = 8, height = 4)
farima(log(ts_export[, 'ex_hommali']), d=1, h=24) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Hommali Rice')
ggsave('../figures/plot_arima_exp_hml.png', width = 8, height = 4)

farima(log(ts_export[, 'ex_pathum']), d=1, h=24) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Pathumtani Rice')
ggsave('../figures/plot_arima_exp_pat.png', width = 8, height = 4)

farima(log(ts_export[, 'ex_white']), d=1, h=24) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'White Rice')
ggsave('../figures/plot_arima_exp_wht.png', width = 8, height = 4)

farima(log(ts_export[, 'ex_glutinous']), d=1, h=24) %>% autoplot() +
  theme_ipsum_rc(base_size = 14) + 
  labs(subtitle = 'Glutinous Rice')
ggsave('../figures/plot_arima_exp_glu.png', width = 8, height = 4)
# Compute CV errors for ETS as e1
e1 <- tsCV(log(data$dom_hommali), fets, h=2, initial = 10)
# Compute CV errors for ARIMA as e2
e2 <- tsCV(log(data$dom_hommali), farima, d=1, h=2, initial = 10)

# MSE
# sqrt(mean(e1[10:13,]^2))
# sqrt(mean(e2[10:13,]^2))

# # MAE
# mean(abs(e1[10:13,]))
# mean(abs(e2[10:13,]))
forecast(ets(log(data$ex_hommali)), h = 2)
forecast(auto.arima(log(data$dom_glutinous), d = 1), h = 2) %>% autoplot()<jupyter_output><empty_output><jupyter_text>VAR<jupyter_code>VARselect(ts_export)
ca.jo(ts_export, ecdet = 'const', K = 4, season = 12) %>% summary()
vecm_export = VAR(diff(ts_export), p = 4, type = 'const', season = 12)
serial.test(vecm_export, lags.pt=12, type="PT.asymptotic")
forecast(vecm_export, h = 24) %>% autoplot()
ggsave('../figures/plot_vecm_exp.png', width = 14, height = 8)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1773482536453035, 'fraction_numerical': 0.03577483892845032, 'mean_word_length': 3.3943006146395978, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 144, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '245420', 'license': 'MIT', 'n_tokens_mistral': 10733, 'n_tokens_neox': 9910, 'n_words': 1824, 'path': 'notebook/backup/demand_forecast-Copy1.ipynb', 'repo_name': 'piyayut-ch/demand_forecast', 'size': '23592'}