n <- 100 # 全時間数
x <- 1:n # 時刻
t <- 2 * x # トレンド（係数：2）
e <- rnorm(n, mean = 0, sd = 10) # 不規則変動
y <- t + e # 合成変動
matplot(x, y, type = 'o', pch = 16, col = 2)
n.tr <- n * 0.8 # 訓練時間（80%）
n.te <- n - n.tr # テスト時間（20%）
ii.tr <- 1:n.tr
ii.te <- (n.tr + 1):n
fit <- lm(y ~ x, data = data.frame(x = x[ii.tr], y = y[ii.tr]))
#summary(fit)
pred <- predict(fit, newdata = data.frame(x = x[ii.te]), interval = 'prediction')
myplot <- function(x, y, pred)
{
matplot(x, y, type = 'o', pch = 16, col = 2, cex = 0.5)
abline(v = n.tr+1, col = 3, lty = 2)
matlines(x[ii.te], pred[, 'fit'], type = 'o', pch = 1, lty = 2, col = 4)
matlines(x[ii.te], pred[, 'lwr'], lty = 3, col = 'lightgray')
matlines(x[ii.te], 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%予測区間'))
}
myplot(x, y, pred)
t <- 0.05 * x^2 - 2 * x # 2次のトレンド
y <- t + e # 合成変動
fit <- lm(y ~ poly(x, 2, raw = T), data = data.frame(x = x[ii.tr], y = y[ii.tr]))
#summary(fit)
pred <- predict(fit, newdata = data.frame(x = x[ii.te]), interval = 'prediction')
myplot(x, y, pred)
T <- 10 # 周期
c <- 30 * sin(2 * pi * x / T) # 周期変動（振幅：30，周期：10）
t <- 2 * x # トレンド（係数：2）
y <- t + c + e # 合成変動
fit <- lm(y ~ poly(x, 2, raw = T), data = data.frame(x = x[ii.tr], y = y[ii.tr]))
#summary(fit)
pred <- predict(fit, newdata = data.frame(x = x[ii.te]), interval = 'prediction')
myplot(x, y, pred)
dir.create('./ts_reg_periodic', showWarnings = F)
yhat <- rep(NA, n)
for (j in 1:T)
{
# j <- 1
ii <- seq(j, n.tr, T)
ii.pred <- ii[length(ii)] + seq(T, T*2, T)
fit <- lm(y ~ x, data = data.frame(x = x[ii], y = y[ii]))
pred <- predict(fit, newdata = data.frame(x = x[ii.pred]), interval = 'prediction')
yhat[ii.pred] <- pred[, 'fit']
png(paste0('./ts_reg_periodic/ts_reg_periodic', j, '.png'))
matplot(x, y, type = 'o', pch = 16, col = 2, cex = 0.5, ylim = c(0, 250))
abline(v = n.tr+1, col = 3, lty = 2)
matlines(x[ii], y[ii], type = 'o', pch = 1, col = 1, cex = 1.5)
matlines(x, yhat, type = 'o', pch = 1, col = 4)
matpoints(x[ii.pred], yhat[ii.pred], pch = 1, col = 1)
dev.off()
#Sys.sleep(2)
}
#system('bash -c 'convert -delay 100 -loop 0 ./ts_reg_periodic/ts_reg_periodic{1..10}.png ./ts_reg_periodic/ts_reg_periodic.gif'')
d0 <- read.csv('https://stats.dip.jp/01_ds/data/load/load_amedas_y2017-2022_tokyo.csv')
head(d0)
## date datetime dow name event md gr year month
## 1 2017-01-01 2017-01-01 00:00:00 日 祝日：元日 年末年始 01-01 6 2017 1
## 2 2017-01-01 2017-01-01 01:00:00 日 祝日：元日 年末年始 01-01 6 2017 1
## 3 2017-01-01 2017-01-01 02:00:00 日 祝日：元日 年末年始 01-01 6 2017 1
## 4 2017-01-01 2017-01-01 03:00:00 日 祝日：元日 年末年始 01-01 6 2017 1
## 5 2017-01-01 2017-01-01 04:00:00 日 祝日：元日 年末年始 01-01 6 2017 1
## 6 2017-01-01 2017-01-01 05:00:00 日 祝日：元日 年末年始 01-01 6 2017 1
## day su mo tu we th fr sa wd ho bw af sp ab mw time hr mi hm rm ws
## 1 1 1 0 0 0 0 0 0 0 1 0 0 1 0 27830 00:00:00 0 0 77.43 0 2.32
## 2 1 1 0 0 0 0 0 0 0 1 0 0 1 0 26340 01:00:00 1 0 71.00 0 2.52
## 3 1 1 0 0 0 0 0 0 0 1 0 0 1 0 25200 02:00:00 2 0 68.86 0 1.52
## 4 1 1 0 0 0 0 0 0 0 1 0 0 1 0 24380 03:00:00 3 0 71.43 0 1.40
## 5 1 1 0 0 0 0 0 0 0 1 0 0 1 0 23890 04:00:00 4 0 74.86 0 1.97
## 6 1 1 0 0 0 0 0 0 0 1 0 0 1 0 23940 05:00:00 5 0 73.71 0 2.25
## tp tp3h tp5h tp3d tp7d
## 1 4.87 5.31 5.53 5.50 6.74
## 2 4.49 5.01 5.20 5.53 6.73
## 3 4.07 4.48 4.90 5.54 6.72
## 4 3.91 4.16 4.60 5.56 6.71
## 5 3.19 3.72 4.11 5.58 6.69
## 6 3.14 3.41 3.76 5.61 6.68
d0$px <- as.POSIXct(d0$datetime)
d0$gw <- d0$mw / 1000 # 測定値 (単位変換：MW -> GW)
d0$gwfit <- NA
d0$gwlwr <- NA
d0$gwupr <- NA
# 訓練データの取得期間フラグ（ここでは予測日の前後１ヶ月間をとった）
is.67 <- '06-01' <= d0$md & d0$md <= '07-31'
# 訓練期間フラグ
is.tr <- '2017-06-01' <= d0$date & d0$date <= '2022-06-30'
# テスト期間フラグ
is.te <- '2022-07-01' <= d0$date & d0$date <= '2022-07-03'
# グラフ表示期間フラグ
is.plot <- '2022-06-28' <= d0$date & d0$date <= '2022-07-03'
d.plot <- d0[is.plot, ] # グラフ表用データ
matplot(d0$px, d0$gw,
main = '電力需要',
xlab = '年', ylab = '電力需要［GW］',
type = 'l', col = 2, ylim = c(20, 70))
for (hr in 0:23)
{
# 予測時間
#hr <- 10
is.hr <- d0$hr == hr
# 訓練データ
d.tr <- d0[is.hr & is.67 & is.tr, ]
# テストデータ（予測区間の正解データとして精度評価に利用）
d.te <- d0[is.hr & is.te, ]
# フィッティング
# （ここでは，トレンドがなく気温の2次多項式と曜日ダミーのモデルとした）
fit <- lm(gw ~ tp + I(tp^2) + dow, data = d.tr)
summary(fit)
# 予測（予測区間付き）
pred <- predict(fit, newdata = d.te, interval = 'prediction')
# 毎時の予測値を保存
d0$gwfit[d0$px %in% d.te$px] <- pred[, 'fit']
d0$gwlwr[d0$px %in% d.te$px] <- pred[, 'lwr']
d0$gwupr[d0$px %in% d.te$px] <- pred[, 'upr']
# 時間別グラフ
dir.create('./load', showWarnings = F)
png(paste0('./load/load_at', hr, '.png'))
matplot(d0$px[is.hr & is.plot], d0$gw[is.hr & is.plot],
main = paste0(hr, '時の電力需要'),
type = 'n', pch = 16, col = 2, xaxt = 'n',
xlim = range(d0$px[is.plot]), ylim = c(20, 70))
px.ax <- seq(d0$px[1], d0$px[nrow(d0)], by = 60*60*24)
abline(v = px.ax, col = 'lightgray', lty = 2)
abline(v = as.POSIXct('2022-07-01'), col = 3, lty = 2) # 予測開始線
matlines(d0$px[is.plot], d0$gw[is.plot], type = 'o', pch = 1, col = 'lightgray')
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')
matlines(d0$px[is.hr & is.plot], d0$gw[is.hr & is.plot], type = 'o', pch = 16, col = 2)
# 凡例
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$gw)/d.te$gw))*100))
dev.off()
}
#system('bash -c 'convert -delay 100 -loop 0 ./load/load_at{0..23}.png ./load/load_at.gif'')
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 12:00:00 2014-01-01 13:00:00 998.1 1002.4 0 37 37.67 5.3 5.67
## 2 2014-01-01 14:00:00 2014-01-01 15:00:00 998.2 1002.5 0 35 35.00 4.8 4.55
## 3 2014-01-01 08:00:00 2014-01-01 09:00:00 1001.9 1006.3 0 44 45.33 3.0 3.20
## 4 2014-01-01 15:00:00 2014-01-01 16:00:00 998.8 1003.1 0 36 35.50 3.9 3.88
## 5 2014-01-01 09:00:00 2014-01-01 10:00:00 1001.7 1006.1 0 41 42.50 4.0 3.78
## 6 2014-01-01 07:00:00 2014-01-01 08:00:00 1001.8 1006.2 0 46 47.00 2.6 2.90
## ws_max rm rm1h tp tp1h tp5h tp3d tp7d wd.y wd_max missing
## 1 9.3 10 60.00 14.8 14.57 11.95 6.59 6.37 南南西 南南東
## 2 8.9 10 58.02 14.8 15.05 14.12 6.76 6.43 南南西 南南西
## 3 5.5 10 60.00 9.4 8.77 7.27 6.25 6.23 南南西 南
## 4 6.3 10 55.98 14.3 14.53 14.60 6.85 6.47 南 南
## 5 9.1 10 60.00 11.3 10.45 8.17 6.34 6.26 南 南
## 6 5.2 10 33.00 8.0 7.50 6.41 6.18 6.20 南南西 南南西
#View(d)
# グラフ
MAIN <- 'アイスクリーム月間世帯支出'
matplot(d$px, d$y, type = 'o', pch = 16, lty =1, col = 2,
main = MAIN, xlab = '年', ylab = '円')
