d0 <- read.csv('https://stats.dip.jp/01_ds/data/e-stat_expenditure.csv', skip = 11)
# タイムスタンプと支出金額のカラムを探しデータを抽出
j.time <- grep('時間軸.月次..コード', colnames(d0))
j.ice <- grep('アイスクリーム', colnames(d0))
d1 <- d0[, c(j.time, j.ice)]
# 時間コード（10桁）からPOSIX準拠時間を作成
yyyy <- substr(paste(d1[, 1]), 1, 4)
mm <- substr(paste(d1[, 1]), 7, 8)
px <- as.POSIXlt(paste(yyyy, mm, '01 12:00:00'), format = '%Y%m%d %H:%M:%S')
d0 <- data.frame(px, date = format(px, '%Y-%m-%d'), y = d1[, 2])
d.ca <- read.csv('https://stats.dip.jp/01_ds/data/calendar/calendar.csv')
d.amedas <- read.csv('https://stats.dip.jp/01_ds/data/amedas/amedas_Tokyo_1h.csv')
d.amedas$date <- format(as.POSIXct(d.amedas$HS), "%Y-%m-%d") # カラムdateを作成
# 暦，気象とマージ
d0 |> merge(d.ca, by = "date") |> merge(d.amedas, by = "date") -> d
head(d)
## date px y dow name event md
## 1 2014-01-01 2014-01-01 12:00:00 396 水 祝日（元日） 年末年始（01-01） 01-01
## 2 2014-01-01 2014-01-01 12:00:00 396 水 祝日（元日） 年末年始（01-01） 01-01
## 3 2014-01-01 2014-01-01 12:00:00 396 水 祝日（元日） 年末年始（01-01） 01-01
## 4 2014-01-01 2014-01-01 12:00:00 396 水 祝日（元日） 年末年始（01-01） 01-01
## 5 2014-01-01 2014-01-01 12:00:00 396 水 祝日（元日） 年末年始（01-01） 01-01
## 6 2014-01-01 2014-01-01 12:00:00 396 水 祝日（元日） 年末年始（01-01） 01-01
## gr year month day wk su mo tu we th fr sa wd.x ho bw af sp ab
## 1 6 2014 1 1 1 0 0 0 1 0 0 0 0 1 0 0 1 0
## 2 6 2014 1 1 1 0 0 0 1 0 0 0 0 1 0 0 1 0
## 3 6 2014 1 1 1 0 0 0 1 0 0 0 0 1 0 0 1 0
## 4 6 2014 1 1 1 0 0 0 1 0 0 0 0 1 0 0 1 0
## 5 6 2014 1 1 1 0 0 0 1 0 0 0 0 1 0 0 1 0
## 6 6 2014 1 1 1 0 0 0 1 0 0 0 0 1 0 0 1 0
## HS HE hpa hpa0 pr hm hm1h ws ws1h
## 1 2014-01-01 13:00:00 2014-01-01 14:00:00 997.8 1002.1 0 36 36.50 3.9 3.97
## 2 2014-01-01 23:00:00 2014-01-02 00:00:00 1002.5 1006.9 0 47 46.00 1.0 1.35
## 3 2014-01-01 15:00:00 2014-01-01 16:00:00 998.8 1003.1 0 36 35.50 3.9 3.88
## 4 2014-01-01 21:00:00 2014-01-01 22:00:00 1002.1 1006.6 0 58 54.83 0.6 0.73
## 5 2014-01-01 00:00:00 2014-01-01 01:00:00 1003.5 1008.0 0 68 67.17 1.0 1.83
## 6 2014-01-01 05:00:00 2014-01-01 06:00:00 1001.9 1006.3 0 48 50.33 3.4 3.00
## ws_max rm rm1h tp tp1h tp5h tp3d tp7d wd.y wd_max missing
## 1 8.1 10 60.00 15.4 15.02 13.20 6.68 6.40 南南西 南南西
## 2 2.3 0 0.00 8.0 8.43 8.09 7.29 6.54 北 北東
## 3 6.3 10 55.98 14.3 14.53 14.60 6.85 6.47 南 南
## 4 1.0 0 0.00 6.4 7.17 9.89 7.16 6.55 北西 西
## 5 1.9 0 0.00 5.3 4.82 6.27 5.81 6.13 北北西 北西
## 6 4.8 0 0.00 7.2 7.03 5.20 6.03 6.16 南 南南西
tail(d)
## date px y dow name event md gr year month day
## 2563 2022-11-01 2022-11-01 12:00:00 660 火 平日 11-01 3 2022 11 1
## 2564 2022-11-01 2022-11-01 12:00:00 660 火 平日 11-01 3 2022 11 1
## 2565 2022-11-01 2022-11-01 12:00:00 660 火 平日 11-01 3 2022 11 1
## 2566 2022-11-01 2022-11-01 12:00:00 660 火 平日 11-01 3 2022 11 1
## 2567 2022-11-01 2022-11-01 12:00:00 660 火 平日 11-01 3 2022 11 1
## 2568 2022-11-01 2022-11-01 12:00:00 660 火 平日 11-01 3 2022 11 1
## wk su mo tu we th fr sa wd.x ho bw af sp ab HS
## 2563 45 0 0 1 0 0 0 0 1 0 0 0 0 0 2022-11-01 23:00:00
## 2564 45 0 0 1 0 0 0 0 1 0 0 0 0 0 2022-11-01 11:00:00
## 2565 45 0 0 1 0 0 0 0 1 0 0 0 0 0 2022-11-01 12:00:00
## 2566 45 0 0 1 0 0 0 0 1 0 0 0 0 0 2022-11-01 13:00:00
## 2567 45 0 0 1 0 0 0 0 1 0 0 0 0 0 2022-11-01 00:00:00
## 2568 45 0 0 1 0 0 0 0 1 0 0 0 0 0 2022-11-01 01:00:00
## HE hpa hpa0 pr hm hm1h ws ws1h ws_max rm rm1h tp
## 2563 2022-11-02 00:00:00 1015.9 1018.8 0 84 81.33 0.4 1.68 1.9 0 0 13.9
## 2564 2022-11-01 12:00:00 1019.7 1022.6 0 60 60.17 2.1 2.20 4.4 0 0 16.6
## 2565 2022-11-01 13:00:00 1018.5 1021.4 0 60 61.50 2.4 2.98 4.4 0 0 17.1
## 2566 2022-11-01 14:00:00 1018.0 1020.9 0 62 62.17 1.7 2.02 4.0 0 0 17.3
## 2567 2022-11-01 01:00:00 1023.1 1026.0 0 74 76.83 2.3 2.72 3.5 0 0 12.8
## 2568 2022-11-01 02:00:00 1023.2 1026.2 0 74 74.83 2.2 2.10 4.7 0 0 12.2
## tp1h tp5h tp3d tp7d wd.y wd_max missing
## 2563 14.32 14.83 14.41 14.15 東南東 北北西
## 2564 16.55 14.85 14.60 13.96 北 北東
## 2565 16.83 15.60 14.54 13.98 北北西 北北西
## 2566 17.33 16.32 14.49 14.00 北 北
## 2567 12.73 13.51 14.75 13.88 北西 北西
## 2568 12.37 13.19 14.76 13.88 北北西 北
#View(d)
# グラフ
MAIN <- 'アイスクリーム月間世帯支出'
matplot(d$px, d$y, type = 'o', pch = 16, lty =1, col = 2,
main = MAIN, xlab = '年', ylab = '円')
d$gwfit <- NA
d$gwlwr <- NA
d$gwupr <- NA
# 訓練データの取得期間フラグ（ここでは予測日の前後１ヶ月間をとった）
is.67 <- '06-01' <= d$md & d$md <= '07-31'
# 訓練期間フラグ
is.tr <- '2014-01-01' <= d$date & d$date <= '2021-12-31'
# テスト期間フラグ
is.te <- '2022-01-01' <= d$date & d$date <= '2022-11-01'
# グラフ表示期間フラグ
is.plot <- '2022-06-28' <= d$date & d$date <= '2022-11-01'
d.plot <- d[is.plot, ] # グラフ表用データ
d$yfit <- NA
d$ylwr <- NA
d$yupr <- NA
is.tr <- '2014-06-01' <= d$date & d$date <= '2021-12-31'
is.te <- '2022-01-01' <= d$date & d$date <= '2022-11-30'
is.plot <- '2014-06-01' <= d$date & d$date <= '2022-11-30'
d.plot <- d0[is.plot, ]
for (mon in 1:11)
{
is.mon <- d$month == mon
d.tr <- d[is.mon & is.tr, ]
d.te <- d[is.mon & is.te, ]
fit <- lm(y ~ year, data = d.tr)
summary(fit)
pred <- predict(fit, newdata = d.te, interval = 'prediction')
d$yfit[d.te$px %in% d.te$px] <- pred[, 'fit']
d$ylwr[d.te$px %in% d.te$px] <- pred[, 'lwr']
d$yupr[d.te$px %in% d.te$px] <- pred[, 'upr']
matplot(d$px[is.mon & is.plot], d$y[is.mon & is.plot],
main = paste(MAIN, ' (',mon, '月) '), xlab = '年', ylab = '円',
type = 'o', pch = 16, col = 2, ylim = c(300, 2000))
abline(v = as.POSIXct('2022-01-01'), col = 3, lty = 2)
matlines(d.te$px, pred[, 'fit'], type = 'o', pch = 1, lty = 2, col = 4)
matlines(d.te$px, pred[, 'lwr'], lty = 3, col = 'lightgray')
matlines(d.te$px, pred[, 'upr'], lty = 3, col = 'lightgray')
legend('topleft', bg = 0,
lty = c(1, 2, 3), pch = c(16, 1, NA), col = c(2, 4, 'lightgray'),
legend = c('原系列', '予測値', '95%予測区間'))
legend('topright', bg = 0,
legend = sprintf('MAPE: %2.1f%%',
mean(abs((pred[, 'fit'] - d.te$y)/d.te$y))*100))
}
