Un article assez simple aujourd’hui. Je voulais examiner les chiffres mensuels du tourisme aux Samoa. En fait, j’ai commencé à faire cela pour les îles du Pacifique en général, mais les problèmes de gestion des données étaient suffisants pour que je n’aille que jusqu’aux Samoa pour l’instant. Ces résultats suscitent actuellement de l’intérêt car ils constitueraient un indicateur relativement opportun des éventuels dommages économiques causés par la crise du carburant liée à la guerre en Iran. Le nombre de visiteurs, l’inflation et le commerce des marchandises font partie du nombre très restreint de statistiques économiques mensuelles publiées dans cette partie du monde (pour un sous-ensemble sélectionné de pays).

Il se passe deux choses dans ce post :

• obtenir le nombre d’arrivées de visiteurs aux Samoa ; et
• en les ajustant de façon saisonnière.

Ce dernier m’intéresse plus que le premier, mais malheureusement, le premier a pris la plupart du temps.

Gestion des données

Voici à quoi ressemble le nombre de visiteurs, une fois que nous les avons tous regroupés dans un seul bloc de données :

Le Bureau des statistiques de Samoa dispose d’un joli classeur Excel jusqu’en mai 2023, mais à partir de cette date, les données ne sont disponibles que dans la mesure où je peux le voir dans les rapports PDF. Heureusement, ils sont tous disponibles sous forme de liens à partir d’une seule page, mais il semble y avoir soit un problème de compétence de ma part, soit une sorte de blocage sur le site Web qui arrête tout téléchargement systématique de tous, j’ai donc dû télécharger tous les PDF à la main, un par un.

Ensuite Claude m’a aidé à écrire un analyseur pour trouver le numéro d’arrivée des visiteurs dans chaque PDF. En fait, Claude n’était pas vraiment doué pour trouver le bon numéro, mais cela m’a donné un modèle que je pouvais adapter, un peu comme autrefois j’aurais utilisé Stack Overflow. Il y a beaucoup de chiffres dans chaque PDF et nous avons besoin du bon : les arrivées de visiteurs, pas le nombre total d’arrivées (oui, c’est un tiret cadratin, mais je les écris tous à la main). Dans ce cas, il s’avère que l’astuce est que les fichiers PDF incluent tous la phrase « Le nombre total de visiteurs pour le mois sous revue s’élevait à [number]», et heureusement, il n’y a aucune autre utilisation du mot « se tenait » dans le document.

L’autre chose délicate consistait à extraire la date réelle à laquelle chaque PDF faisait référence. Ensuite, il fallait tout tester. Au final, il aurait été plus rapide de saisir manuellement les 35 chiffres dont j’avais besoin. Mais voici le code qui gère toutes les données.

library(tidyverse)
library(rvest)
library(httr2)
library(lubridate)
library(pdftools)   
library(readxl)
library(seasonal)
library(tsibble)
library(fable)
library(feasts)
library(ggtext)
library(scales)

#' Extract date from messy filenames
extract_date <- function(messy_date){
  tibble(path = messy_date) |>
  mutate(
    stem      = tools::file_path_sans_ext(basename(path)),
    # Remove any leading word followed by _ or - (e.g. "Migration_")
    stem      = str_remove(stem, "^([A-Za-z]+[_-])+(?=[A-Z])"),
    month_str = str_extract(stem, "[A-Za-z]{3,}"),
    # Normalise non-standard abbreviations
    month_str = str_replace(month_str, "^Sept$", "Sep"),
    year_str  = str_extract(stem, "\\d{4}|\\d{2}"),
    year      = if_else(nchar(year_str) == 2,
                        as.integer(paste0("20", year_str)),
                        as.integer(year_str)),
    date      = as.Date(paste(year, month_str, "01"), format = "%Y %B %d") |>
                  coalesce(as.Date(paste(year, month_str, "01"), format = "%Y %b %d"))
  ) |>
  pull(date)
}

test <- c("samoa_pdfs/April_25.pdf", "samoa_pdfs/Feb_25.pdf", "samoa_pdfs/Feb_26.pdf", 
"samoa_pdfs/Jan_26.pdf", "samoa_pdfs/January_25.pdf", "samoa_pdfs/July_25.pdf", 
"samoa_pdfs/June-25.pdf", "samoa_pdfs/March_2025.pdf", "samoa_pdfs/March_2026.pdf", 
"samoa_pdfs/May_25.pdf", "samoa_pdfs/Migration_April-2026.pdf",
"samoa_pdfs/Sept-24.pdf", "samoa_pdfs/Migration_Rep_June_2023.pdf")

extract_date(test)

#-------------------PDFs for recent data------------------
# For the more recent years no Excel tables are published, so need
# to use the PDFs and extract total from there
# These had to be downloaded by hand - nothing I tried was able to automate
# that. Download from  and save
# in a subfolder /samoa_pdfs/.

pdf_dir <- "samoa_pdfs"
tbl <- tibble(local_path = list.files(pdf_dir, pattern = ".pdf$", full.names = TRUE))

parse_pdf_visitors <- function(path) {
  txt <- tryCatch(pdf_text(path), error = function(e) {
    message("  [WARN] pdftools could not read: ", basename(path))
    NULL
  })
  if (is.null(txt)) return(NA_integer_)
  full_text <- paste(txt, collapse = "\n")
  patterns <- c(
    "stood at [^0-9(]{0,10}\\(?([0-9,]+)\\)?"
  )
  for (pat in patterns) {
    m <- str_match(full_text, pat)[, 2]
    if (!is.na(m)) return(as.integer(str_remove_all(m, ",")))
  }
  message("  [WARN] Could not parse visitor count from: ", basename(path))
  NA_integer_
}

found <- tbl |> 
  filter(file.exists(local_path))

message("  Parsing ", nrow(found), " local PDFs...")
pdf_tbl <- found |>
  mutate(visitors = map_int(local_path, \(p) {
    message("  ", basename(p))
    parse_pdf_visitors(p)
  })) |>
  filter(!is.na(visitors)) |>
  mutate(date = extract_date(local_path))

#-----------------Excel versions for older data------------
# For May 2023 and earlier we can get the data for multiple months
# at a time from Table 1 of the Excel tables. The May 2023 Excel
# file goes back to 2017 January (although the rows are hidden)

fn <- "May_23.xlsx"
if(!file.exists(fn)){
  download.file("
                destfile = fn, mode = "wb")
}
x <- read_excel(fn, sheet = "Table 1", 
                 range = "D48:D130",
                 col_names = "visitors") |> 
  drop_na() |> 
  pull(visitors)

historical <- tibble(visitors = x, 
               date = seq(as.Date("2017-01-01"), as.Date("2023-05-01"), by = "month"))

#-----------combine and test-----------------------
samoa_visitors <- pdf_tbl |> 
  select(date, visitors) |>
  bind_rows(historical) |> 
  arrange(date) |> 
  mutate(date_month = yearmonth(date)) |> 
  as_tsibble(index = date_month)


# Test - some hand picked test cases, 4 from PDFs and 3 from the Excel
samoa_test <- tribble(~date, ~correct_visitors,
                      "2023-08-01", 16471,
                      "2024-04-01", 12644,
                      "2024-08-01", 17248,
                      "2026-04-01", 14188,
                      "2018-02-01", 7413,
                      "2020-12-01", 195,
                      "2022-06-01", 866) |> 
  mutate(date = as.Date(date))

stopifnot(
  samoa_test |> 
    anti_join(samoa_visitors, by = c("date", "correct_visitors" = "visitors")) |> 
    nrow() == 0
)

Et voici le code qui dessine le graphique de série chronologique de base que j’ai utilisé plus tôt :

the_caption = "Source: Samoa Bureau of Statistics"

ggplot(samoa_visitors, aes(x = date, y = visitors)) +
  geom_line() +
  scale_y_continuous(label = comma) +
  labs(x = "",
       y = "Visitor arrivals",
       title = "Visitor arrivals per month to Samoa",
       subtitle = "Unadjusted originals",
       caption = the_caption)

Modélisation de la décomposition saisonnière

Ouf. Ok, passons à la partie la plus amusante du mannequinat. J’ai l’intention d’utiliser l’outil X-13ARIMA-SEATS. X‑13ARIMA‑SEATS est un programme développé et maintenu par le US Census Bureau pour la désaisonnalisation et la décomposition des séries chronologiques. Il s’adaptera automatiquement à un modèle de série chronologique SARIMA (moyenne mobile intégrée autorégressive saisonnière), intègre des méthodes d’identification et de traitement des valeurs aberrantes, s’ajuste par défaut au nombre de jours de bourse et aux vacances mobiles de Pâques, et permet à l’utilisateur de spécifier des variables explicatives de régression supplémentaires si vous le souhaitez. C’est la référence mondiale en matière de désaisonnalisation des statistiques officielles.

X-13ARIMA-SEATS est disponible en R via le seasonal paquet (par Christoph Sax et Dirk Eddelbuettel). En aval de cela, le feasts et fable Les packages de (Mitchell O’Hara-Wild, Rob Hyndman et Earo Wangmake) facilitent le travail dans une approche tabulaire et ordonnée.

La période Covid est un facteur dominant évident dans le tourisme au cours des dernières décennies, et cela apparaît dans le premier graphique que j’ai présenté ci-dessus. Je m’intéresse également à la période écoulée depuis le début de la guerre entre les États-Unis, Israël et l’Iran, pour voir si cela a un impact. Il n’y a que deux mois de données (mars et avril 2026) depuis le début de la guerre, il faudrait donc qu’un impact soit dramatique pour se manifester, mais cela vaut la peine de vérifier.

X-13ARIMA-SEATS ajustera par défaut les prévisions lorsque vous lui demanderez de modéliser. Vous devrez donc fournir les valeurs des variables du régresseur x pour couvrir non seulement la période des données, mais aussi quelques périodes à venir. Je fais cela sous forme de simples vecteurs de séries chronologiques. Je n’ai pas encore trouvé le moyen simple de le faire dans le monde tabulaire de fable. Heureusement, il semble que je puisse utiliser ces vecteurs plus tard, même dans fable. Dans un instant, j’utiliserai les deux seasonal directement et via fable pour être sûr d’obtenir les mêmes résultats. En fait, mettre à jour mes connaissances datées de la modélisation de séries chronologiques pour utiliser fable était l’une de mes principales motivations pour tout cet exercice.

Voici le code qui crée les x régresseurs que j’utiliserai pour la présence de la pandémie de Covid et pour la guerre en Iran :

#----------------x regressor variables-----------------

# Covid time series indicator to use as a regressor
covid_reg <- ts(
  as.numeric(seq(as.Date("2017-01-01"), as.Date("2029-04-01"), by = "month") %in%
    seq(as.Date("2020-04-01"), as.Date("2022-07-01"), by = "month")),
  start     = c(2017, 1),
  frequency = 12
)

# Iran war time series indicator
war_reg <- ts(
  as.numeric(seq(as.Date("2017-01-01"), as.Date("2029-04-01"), by = "month") %in%
    seq(as.Date("2026-03-01"), as.Date("2026-07-01"), by = "month")),
  start     = c(2017, 1),
  frequency = 12
)

Directement avec seasonal

Ok, c’est l’heure du mannequinat. Premièrement, en utilisant simplement un vecteur de séries chronologiques à l’ancienne sur les arrivées de visiteurs aux Samoa et les seasonal package directement, voici l’installation du modèle X-13ARIMA-SEATS avec des valeurs par défaut utilisant à la fois les régresseurs de guerre Covid et Iran :

sa_ts <- ts(samoa_visitors$visitors, frequency = 12, start = c(2017, 1))

fit_ts_war <- seas(sa_ts, xreg = cbind(covid_reg, war_reg))
summary(fit_ts_war)

Hyper simple. Cela nous donne ce résultat :

Call:
seas(x = sa_ts, xreg = cbind(covid_reg, war_reg))

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
xreg1             -3.10069    0.14567 -21.286  < 2e-16 ***
xreg2             -0.15937    0.13098  -1.217    0.224    
LS2020.Mar        -0.85036    0.14567  -5.838 5.29e-09 ***
AO2020.Dec        -0.75529    0.12350  -6.115 9.63e-10 ***
LS2021.May        -0.79012    0.13233  -5.971 2.36e-09 ***
AO2021.Jul        -1.36973    0.12378 -11.066  < 2e-16 ***
LS2022.May         0.89068    0.13233   6.731 1.68e-11 ***
LS2022.Aug        -1.07689    0.19608  -5.492 3.97e-08 ***
AR-Nonseasonal-01 -0.59006    0.07041  -8.380  < 2e-16 ***
MA-Seasonal-12     0.99961    0.08044  12.427  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SEATS adj.  ARIMA: (1 1 0)(0 1 1)  Obs.: 112  Transform: log
AICc:  1608, BIC:  1634  QS (no seasonality in final):3.295  
Box-Ljung (no autocorr.): 26.42   Shapiro (normality): 0.9721 *

Quelques éléments clés à noter ici.

Les données ont été transformées en log, ce qui est une bonne chose. J’aurais certainement choisi de faire cela, ou du moins une transformation en racine carrée, étant donné la façon dont la variance des arrivées de visiteurs augmente comme sa moyenne, partout dans le monde.

Le modèle finalement adopté est décrit comme ARIMA (1 1 0)(0 1 1). Cela signifie la série principale comme terme d’autorégression sur un décalage, après un cycle de différenciation ; et la partie saisonnière a un terme de moyenne mobile en retard, après un cycle de différenciation. Il s’agit d’un modèle tout à fait normal pour les chiffres du tourisme. cela indique une tendance/impulsion générale (la différence dans la série principale), une tendance pour les valeurs d’un mois à être liées d’une manière ou d’une autre à celles du mois précédent (dans ce cas avec une corrélation négative de -0,59) et un fort effet de saisonnalité annuelle qui évolue lentement dans le temps.

La variable muette Covid (xreg1) est fortement négativement significatif. Avec un coefficient de -3,1 et la transformée logarithmique de la variable réponse, cela signifie que pendant la période Covid les arrivées réelles étaient en moyenne de exp(-3.1) = 0.045 (soit 4,5% soit en baisse de 95,5%) pendant les périodes hors Covid.

En revanche, nous n’avons pas d’effet statistiquement significatif pour la guerre en Iran (xreg2). Pour une analyse ultérieure, je supprimerai ce régresseur x car je ne veux pas qu’il complique la tendance récente et les chiffres désaisonnalisés.

Les données de six mois ont été identifiées comme valeurs aberrantes et contrôlées de manière appropriée. Tout cela se situe dans la période Covid difficile à modéliser de 2020 à 2022.

Un point à noter est que nous n’avons pas d’effet Pâques. Je suis sûr à 100 % qu’il y a réellement un effet de Pâques dans les arrivées de visiteurs aux Samoa, mais 9 années de données ne suffisent pas pour le montrer. Pâques a parfois lieu en avril et parfois en mars. Mais depuis 2017, cela se produit chaque année en avril, sauf en 2024. Ce n’est tout simplement pas suffisamment de variation pour le distinguer des impacts saisonniers mensuels réguliers.

Pour vérifier si je me souviens que Pâques est bien cochée par défaut dans X-13ARIMA-SEATS, j’adapte un modèle aux données bien connues des compagnies aériennes Box et Jenkins :

> # Another examplefor comparison
+ m <- seas(AirPassengers)
+ summary(m)

Call:
seas(x = AirPassengers)

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
Weekday           -0.0029497  0.0005232  -5.638 1.72e-08 ***
Easter[1]          0.0177674  0.0071580   2.482   0.0131 *  
AO1951.May         0.1001558  0.0204387   4.900 9.57e-07 ***
MA-Nonseasonal-01  0.1156204  0.0858588   1.347   0.1781    
MA-Seasonal-12     0.4973600  0.0774677   6.420 1.36e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SEATS adj.  ARIMA: (0 1 1)(0 1 1)  Obs.: 144  Transform: log
AICc: 947.3, BIC: 963.9  QS (no seasonality in final):    0  
Box-Ljung (no autocorr.): 26.65   Shapiro (normality): 0.9908 

On voit ici que le nombre de jours de la semaine dans un mois et le congé de Pâques décalé sont bien dans le modèle final, sans que j’aie dû demander leur vérification. Pâques a un impact positif sur cette série de passagers aériens (1949 à 1960), et le nombre de jours de la semaine dans un mois a un impact négatif moindre.

Ainsi, après cette digression de vérification de Pâques, je réajuste le modèle pour mes arrivées de visiteurs aux Samoa sans le régresseur de guerre et j’obtiens un résultat essentiellement identique :

> fit_ts <- seas(sa_ts, xreg = covid_reg)
> summary(fit_ts)

Call:
seas(x = sa_ts, xreg = covid_reg)

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
xreg              -3.10045    0.14690 -21.106  < 2e-16 ***
LS2020.Mar        -0.83292    0.14617  -5.698 1.21e-08 ***
AO2020.Dec        -0.75463    0.12449  -6.062 1.35e-09 ***
LS2021.May        -0.78975    0.13356  -5.913 3.36e-09 ***
AO2021.Jul        -1.36737    0.12476 -10.960  < 2e-16 ***
LS2022.May         0.89113    0.13356   6.672 2.52e-11 ***
LS2022.Aug        -1.07733    0.19782  -5.446 5.15e-08 ***
AR-Nonseasonal-01 -0.58577    0.07071  -8.284  < 2e-16 ***
MA-Seasonal-12     0.99912    0.07827  12.765  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SEATS adj.  ARIMA: (1 1 0)(0 1 1)  Obs.: 112  Transform: log
AICc:  1607, BIC:  1631  QS (no seasonality in final):2.501  
Box-Ljung (no autocorr.): 26.63   Shapiro (normality): 0.9706 *

Avec fable et feasts

tsibble, fable et feasts sont un brillant ensemble de packages qui vous permettent de travailler avec des séries chronologiques dans R de manière plus tabulaire et tidyversede manière conviviale que les différentes structures de données de séries chronologiques plus anciennes ne vous le permettent. Cependant, la plupart de mon travail avec les séries chronologiques a eu lieu avant qu’elles ne soient disponibles, donc je manque de confiance dans leur fonctionnement. Heureusement, cela semble assez simple.

J’avais déjà franchi l’étape critique plus tôt avec as_tsibble(index = date_month)en précisant mon principal samoa_visitors tibble est en fait un tibble de série chronologique. Maintenant, la modélisation utilisant ce tsibble est assez simple :

#----------fable version-------
fit_fb <- samoa_visitors |> 
  model(X_13ARIMA_SEATS(
    visitors ~ xreg(covid_reg)
  ))

report(fit_fb)

Notez que nous utilisons report() plutôt que summary() pour obtenir le résultat final. Je ne vais pas l’imprimer ici car il est littéralement identique à ce que nous avons obtenu plus tôt. summary(fit_ts).

Le principal attrait pour moi dans l’utilisation du fable/feasts L’approche est qu’elle correspond mieux à la fois à mon flux de travail de gestion des données et à mon approche de ggplot2 graphique. voici donc une belle décomposition de la série temporelle originale, produite avec autoplot(

fit_fb |> 
  components() |> 
  autoplot() 

Il convient de noter que dans cette décomposition, la série initiale des « visiteurs » et la tendance sont sur l’échelle originale, mais les composantes « saisonnières » et « irrégulières » sont exprimées sous forme de multiplicateurs. Ainsi, une valeur saisonnière de 1,2 signifie que, pour un mois donné, la valeur est 20 % plus élevée en raison de la saisonnalité qu’autrement.

Et voici ma version finale de présentation des données :

Produit avec ce code :

comp_data <- fit_fb |> 
  components() |> 
  # the 'trend' that comes straight from the decomposition does
  # not adjust ofr the Covid coefficient and looks pretty weird
  # so it is more intuitive to present it after adjustment for
  # Covid, which we need to calculate by multiplying (because of
  # the log transform that SEATS used autoamtically):
  mutate(covid = as.numeric(covid_reg)[1:nrow(samoa_visitors)],
         trend_adj = trend * exp(covid * coef(fit_ts)[1])) 

comp_data|> 
  ggplot(aes(x = date_month, y = season_adjust)) +
  geom_line(linewidth = 1.3) +
  geom_line(aes(y = trend_adj), colour = "steelblue", alpha = 0.9, linewidth = 1.2) + 
  geom_point(aes(y = visitors), colour = "grey70") +
  scale_y_continuous(label = comma) +
  labs(x = "",
       y ="Visitor arrivals",
       title = "Visitor arrivals per month to Samoa",
       subtitle = "<span style="color:grey60">Original</span>, seasonally adjusted <span style="color:steelblue">and trend (adjusted for Covid period).</span>",
       caption = "Source: data from Samoa Bureau of Statistics. Seasonal adjustment by freerangestats.info.") +
  theme(plot.subtitle = element_markdown())

Après tout cela, quelle perception avons-nous ? En fait, pas grand-chose en réalité, à part les tendances aveuglantes et évidentes de dévastation de l’industrie causées par Covid et la tendance à la croissance lente et bruyante depuis lors. Nous sommes au moins bien placés pour commenter l’impact de la crise du carburant sur le tourisme, et pouvons dire qu’il n’y a pas encore de preuve évidente. Si et quand nous constatons l’impact, nous pourrons en parler en termes de tendances et de variations aléatoires, après avoir supprimé l’élément saisonnier. C’est donc utile.

Eh bien, c’est tout. Peut-être que dans un article ultérieur j’ajouterai les autres pays insulaires du Pacifique avec des données mensuelles sur le tourisme – Fidji, Vanuatu, les Îles Cook et la Polynésie française étant les principaux pays que je connais.

En rapport