観光来訪者数の時系列分析

伊東宏樹

2025-05-10

はじめに

  • 観光客の人数の時系列変化は基本データ

  • 金沢市への観光来訪者数がどのように変化しているのか

    • 全般的な傾向

    • 能登半島地震の影響

    • 季節によるパターン

時系列分析

  • 時系列データは
    • 通常のデータと同様に扱えない場合がおおい
    • 自己相関がある → 独立ではない

何も考えずに時間を説明変数に回帰すると

ランダムウオークでもだいたい有意（p<0.05）になる

状態空間モデル

  • 時系列データをあつかうモデルのひとつ

  • 観測値は、直接観測することのできない「状態」から生成されると考える

  • 状態は、時系列に沿って変化する

グラフィカルモデル

x[t]: 時間tにおける状態, y: 時間tにおける観測値

観光来訪者数データ

  • 公益社団法人日本観光振興協会「デジタル観光統計オープンデータ」の「市区町村観光来訪者数」2021〜2024年を使用
  • スマートフォンの位置情報データを利用して収集
  • 観光来訪者数: 「推定発地から半径20km以上離れた調査地点に滞在した者。但し、調査地点勤務者を除く。」（デジタル観光統計オープンデータの概要）
    • 調査地点とは、調査対象として設定された観光地点

データ読み込み

ダウンロードしたCSVファイルを読み込み、1つのデータフレームにまとめる

> data_dir <- "data"
> data_files <- c("city2021.csv", "city2022.csv",
+                 "city2023.csv", "city2024.csv")
> 
> data <- purrr::map(data_files,
+           \(f) read_csv(file = file.path(data_dir, f),
+                         locale = locale(encoding = "CP932"))) |>
+   dplyr::bind_rows()

データの出典: デジタル観光統計オープンデータ（https://www.nihon-kankou.or.jp/home/jigyou/research/d-toukei/）（2025年3月11日ダウンロード）

データの抽出

金沢市のデータを抽出

> data_Kanazawa <- data |>
+   dplyr::filter(`地域名称` == "金沢市")

データの確認

> print(data_Kanazawa)
# A tibble: 48 × 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 金沢市  
 7  2021     7 市区町村 観光来訪者数……             17 石川県          17201 金沢市  
 8  2021     8 市区町村 観光来訪者数……             17 石川県          17201 金沢市  
 9  2021     9 市区町村 観光来訪者数……             17 石川県          17201 金沢市  
10  2021    10 市区町村 観光来訪者数……             17 石川県          17201 金沢市  
# ℹ 38 more rows
# ℹ 1 more variable: 人数 <dbl>

データの加工

あとの可視化と解析のためデータを加工する

  • 年と月から、“年/月”となる文字列変数YMをつくる
  • 2021年1月からの月数をMonthsに
  • 人数を1万人単位にしてNumに
  • 必要な列（変数）だけ残す
> data_Kanazawa <- data_Kanazawa |>
+   dplyr::mutate(YM = str_c(`年`, `月`, sep = "/"),
+                 Months = (`年` - 2021) * 12 + `月`,
+                 Num = `人数` / 10000,
+                 .keep = "none")

データをグラフに

状態空間モデルの構築

  • システムモデル

    • 状態の変化

  • 観測モデル

    • 状態から観測値が得られる過程

システムモデル

以下の3つの成分からなるモデル

  • レベル（+トレンド）: レベルは、下の回帰成分と季節成分を除いた変化に対応
    • あるトレンド（傾き）をもって変化し、トレンドも変化する
  • 季節成分: 月ごとの変動（12か月周期）
  • 回帰成分: 能登半島地震の影響を検討するため、地震発生の2024年1月以降を示すダミー変数を説明変数としてモデルに含める

それぞれ、あるノイズの値が加わって前の月から変化

観測モデル

  • レベル＋季節＋回帰に、月ごとのノイズ（たまたま多かったり、少なかったり）が加わるとする

今回のデータでは観測値はもともとは計数値だが、数が大きく、1万人単位としているので、観測モデルのノイズは正規分布するものとして扱う

KFASパッケージ

Rで状態空間モデルをあつかうパッケージはいくつかあるが、今回はKFASを使用

  • 状態空間モデルの定義と観測データへのあてはめ

  • カルマンフィルタによるフィルタリング・平滑化・予測

  • 正規分布以外（ポアソン分布など）の観測データにも対応

モデルの定義とあてはめ

SSModel関数でモデルを定義する。引数のQとHはそれぞれシステムモデルと観測モデルのノイズの分散。fitSSM関数であてはめを実行。

> 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 = data_Kanazawa)
> fit <- fitSSM(model, inits = c(0, 0, 0, 0, 0))

平滑化

KFS関数によりカルマンフィルタを適用し、平滑化した状態の値を取得

> smooth <- KFS(fit$model) |>
+   coef(filtered = FALSE)

レベル成分

地震の影響と季節成分を除いた数に対応

観測値と重ねると

トレンド

1か月あたりの増加数

季節成分

能登半島地震の影響

予測

地震の影響はなくなったものとして12か月先までを予測

まとめ

  • 能登半島地震の影響による来訪者数の減少は1月あたり4.4万人程度と推定された
    • 厳密に言うと、2024年1月以降の減少なので、同時期の他の影響もあるかもしれない
  • 2021年の秋以降、来訪者数は増加傾向ですが、増加率は最近減少
  • 来訪者数は毎年11月にもっとも多く、つづいて5月、3月の順
  • 来訪者数が少なくなるのは1月と2月