library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(KFAS)
## Please cite KFAS in publications by using:
##
## Jouni Helske (2017). KFAS: Exponential Family State Space Models in R. Journal of Statistical Software, 78(10), 1-39. doi:10.18637/jss.v078.i10.金沢市中心部を歩いていると、多くの観光客を目にします。観光客の数の変化を時系列分析して、どのような傾向があるのかを検討しました。
使用パッケージ
まずは使用するパッケージを読み込みます。データ処理にはtidyverseを使います。
時系列分析のために今回はKFASパッケージを利用して、状態空間モデルのあてはめをおこないます。
データ
観光客のデータとして、公益社団法人日本観光振興協会「デジタル観光統計オープンデータ」の「市区町村観光来訪者数」2021〜2024年を使用しました。
このデータは、スマートフォンの位置情報データを利用して収集されたもので、「観光来訪者数」は「推定発地から半径20km以上離れた調査地点に滞在した者。但し、調査地点勤務者を除く。」と定義されています（デジタル観光統計オープンデータの概要）。調査地点とは、調査対象として設定された観光地点です。
ダウンロードしたCSVファイル（city2021.csv, city2022.csv, city2023.csv, city2024.csv）は、dataフォルダに入れておきます。これらのファイルを読み込み、読み込んだデータを1つのデータフレームにまとめておきます。
data_dir <- "data"
data_files <- c("city2021.csv", "city2022.csv", "city2023.csv", "city2024.csv")
data_list <- purrr::map(data_files,
\(f) read_csv(file = file.path(data_dir, f),
locale = locale(encoding = "CP932"),
show_col_types = FALSE))
data <- dplyr::bind_rows(data_list)データの出典: デジタル観光統計オープンデータ（https://www.nihon-kankou.or.jp/home/jigyou/research/d-toukei/）（2025年3月11日ダウンロード）
データの抽出と加工
金沢市のデータを抽出します。
k_data <- data |>
dplyr::filter(`地域名称` == "金沢市")以下のようなデータとなっています。
head(k_data)
## # A tibble: 6 × 9
## 年 月 地域区分 データ区分 都道府県コード 都道府県名 地域コード 地域名称
## <dbl> <dbl> <chr> <chr> <dbl> <chr> <dbl> <chr>
## 1 2021 1 市区町村 観光来訪者数…… 17 石川県 17201 金沢市
## 2 2021 2 市区町村 観光来訪者数…… 17 石川県 17201 金沢市
## 3 2021 3 市区町村 観光来訪者数…… 17 石川県 17201 金沢市
## 4 2021 4 市区町村 観光来訪者数…… 17 石川県 17201 金沢市
## 5 2021 5 市区町村 観光来訪者数…… 17 石川県 17201 金沢市
## 6 2021 6 市区町村 観光来訪者数…… 17 石川県 17201 金沢市
## # ℹ 1 more variable: 人数 <dbl>あとの可視化と解析のためデータを加工します。
年と月から、“年/月”となる文字列変数YMをつくります。2021年1月からの月数を
Monthsという変数に格納します。人数を1万人単位にして
Numという変数に格納します。必要な列（変数）だけ残します。
k_data <- k_data |>
dplyr::mutate(YM = str_c(`年`, `月`, sep = "/"),
Months = (`年` - 2021) * 12 + `月`,
Num = `人数` / 10000,
.keep = "none")可視化
グラフで表示します。
max_months <- max(k_data$Months)
p <- ggplot(k_data) +
geom_line(aes(x = Months, y = Num), linewidth = 1) +
scale_x_continuous(name = "年/月",
breaks = seq(1, max_months - 5, 6),
minor_breaks = seq(1, max_months, 3),
labels = k_data$YM[seq(1, max_months - 5, 6)]) +
scale_y_continuous(name = "人数（万人）") +
theme_gray(base_family = "YuGothic")
plot(p)
新型コロナウイルス感染症対策のまん延防止等重点措置が終了したのが2022年3月21日ですが、だいたいそれ以降、観光来訪者数が増加しているように見えます。
状態空間モデル
モデル
KFASパッケージを使用して状態空間モデルを作成し、データにあてはめます。状態空間モデルは、実際には観測できないシステムの状態と、状態から観測値が得られる過程とをそれぞれ明示的にモデル化したものです。
状態モデルは以下の成分から構成されています。
回帰: 能登半島地震の影響を検討するため、地震発生の2024年1月以降を示すダミー変数
quakeを説明変数としてモデルに含めます。地震の影響は2024年1〜12月の間ずっと同じだけあると仮定しています。トレンド: 次数2のモデル（水準成分と傾き成分）とします。
季節成分: 毎月の変動をダミー変数で含めます。
新型コロナウイルスの影響ですが、だんだん影響が薄れると想定し、トレンド成分に入ることを考えて、とくに回帰とはしませんでした。
また、目的変数Numはもともとは計数値ですが、数が大きく、1万人単位としていますので、観測モデルの誤差は正規分布するものとして扱います。
作成したモデルをmodelオブジェクトに入れ、fitSSM関数であてはめます。引数のQとHはそれぞれ状態モデルと観測モデルの分散で、これを推定します。
quake <- c(rep(0, 36), rep(1, 12))
model <- SSModel(Num ~ SSMregression(~ quake, Q = NA) +
SSMtrend(degree = 2,
Q = list(matrix(NA), matrix(NA))) +
SSMseasonal(period = 12, Q = NA, sea.type = "dummy"),
H = NA, data = k_data)
fit <- fitSSM(model, inits = c(0, 0, 0, 0, 0))これくらい複雑なモデルになると、分散の初期値（inits引数）によっては局所最適にはまるようなので、初期値をいくつか変えて対数尤度が最大になるものを選びました。本来ならグリッドサーチとかするべきなのかもしれませんが、そこまではやっていません。
logLik(fit$model)
## [1] -130.5529カルマンフィルタにより平滑化し、その値を取得します。
smooth <- KFS(fit$model) |>
coef(filtered = FALSE)地震の影響
地震影響の係数の推定値をみてみます。
smooth[, "quake"][1]
## [1] -4.377758その標準偏差も確認します。
sqrt(fit$model$Q[1, 1, 1])
## [1] 0.008701566標準偏差からみて、係数は有意に負で、1か月の来訪者数がおよそ4.4万人減少したと推定されました。
各成分
平滑化した水準成分を図示します。地震の影響と季節成分を除いた数に対応します。
ggplot(data = data.frame(Months = 1:max_months,
level = c(smooth[, "level"]))) +
geom_line(mapping = aes(x = Months, y = level),
color = "red", linewidth = 1) +
scale_x_continuous(name = "年/月",
breaks = seq(1, max_months - 5, 6),
minor_breaks = seq(1, max_months, 3),
labels = k_data$YM[seq(1, max_months - 5, 6)]) +
theme_gray(base_family = "YuGothic")
2021年の秋以降急増して、最近は頭打ち傾向のように見えます。
つづいて、傾き成分です。
ggplot(data = data.frame(Months = 1:max_months,
slope = c(smooth[, "slope"]))) +
geom_line(mapping = aes(x = Months, y = slope),
color = "red", linewidth = 1) +
scale_x_continuous(name = "年/月",
breaks = seq(1, max_months - 5, 6),
minor_breaks = seq(1, max_months, 3),
labels = k_data$YM[seq(1, max_months - 5, 6)]) +
theme_gray(base_family = "YuGothic")
2022年以降、傾き（増加率）が減少しているようです。2024年10月以降は0.4くらいで横ばいになっていますが、これが実際にそうなのかはもっと後のデータも含めて検討する必要があるでしょう。
最後に季節成分です。
ggplot(data = data.frame(Months = 1:max_months,
seasonal = c(smooth[, "sea_dummy1"]))) +
geom_line(mapping = aes(x = Months, y = seasonal),
color = "red", linewidth = 1) +
scale_x_continuous(name = "年/月",
breaks = seq(1, max_months - 5, 6),
minor_breaks = seq(1, max_months, 1),
labels = k_data$YM[seq(1, max_months - 5, 6)]) +
theme_gray(base_family = "YuGothic")
毎年11月に来訪者数が多く、つづいて5月、3月となっています。一方、1月と2月の落ち込みが激しくなっています。
まとめ
金沢市の2021〜2024年の観光来訪者数を状態空間モデルで解析した結果、以下のことが推測されました。
- 能登半島地震の影響による来訪者数の減少は1月あたり4.4万人程度（厳密に言うと、2024年1月以降の減少なので、同時期の他の影響もあるかもしれません）
- 2021年の秋以降、来訪者数は増加傾向にあるが、増加率は最近減少している
- 来訪者数は毎年11月にもっとも多く、つづいて5月、3月の順
- 来訪者数が少なくなるのは1月と2月