Analiza.Rmd

---
title: "Modelowanie szeregów czasowych z wykorzystaniem modeli ARIMA"
author: "Krawiec Piotr"
date: "18/06/2021"
output: 
  md_document:
    variant: markdown_github
  pdf_document:
    number_sections: true
    keep_tex: yes
csl: cit.csl
lang: pl
abstract: |
  Celem pracy jest analiza dwóch szeregów czasowych z wykorzystaniem języka R. Pierwszego zawierającego silną sezonowość, drugiego silny trend. Szeregi poddane dekompozycji, usunięciu trendu i sezonowości. Zbadana i potwierdzona została stacjonarność reszt obu szeregów, a ostatecznie zostały modelowane modelem ARIMA. Przewidywania dotyczące dalszego ich przebiegu zostały stworzone z pomocą metod naiwnych uwzględniających dryf, a także sezonowość. Reszty żadnego z zaprezentowanych modeli nie przechodzą testów.
---

```{=tex}
\newpage
\tableofcontents
\newpage
```

------------------------------------------------------------------------

```{r setup, include=FALSE, warning=FALSE}
library(knitr)
hook_output <- knit_hooks$get("output")
knit_hooks$set(output = function(x, options) {
  lines <- options$output.lines
  if (is.null(lines)) {
    return(hook_output(x, options))  # pass to default hook
  }
  x <- unlist(strsplit(x, "\n"))
  more <- "..."
  if (length(lines)==1) {        # first n lines
    if (length(x) > lines) {
      # truncate the output, but add ....
      x <- c(head(x, lines), more)
    }
  } else {
    x <- c(more, x[lines], more)
  }
  # paste these lines together
  x <- paste(c(x, ""), collapse = "\n")
  hook_output(x, options)
})
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(cache = TRUE)
Sys.setlocale("LC_ALL", "Polish")
```

# Szereg - Rozwój biznesu

Na szereg ten składają się dane po chodzące ze strony [FRED](https://fred.stlouisfed.org/series/BUSAPPWNSAUS "FRED Economic data - Business Applications for United States").Dane zbierane są przez U.S Census Bureau, obejmują lata 2006-2021. Zbierane są w tygodniowych odstępach i dotyczą ilości wniosków o wydanie identyfikatora EAN (Employer Identyfication Number). Każdy pracodawca, korporacja, organizacja non-profit itp muszą posiadać takie numery, aby móc rozliczać się z podatku. Jest to zatem dobry wskaźnik tego ile nowych biznesów powstaje.

Do korzyści jakie przyniesie prognoza należy przewidywanie rozwoju gospodarki, gdyż nowo powstające biznesy mogą świadczyć o tym że w kraju panują korzystne warunki do rozwoju biznesu. Analiza szeregu pozwoli też przewidzieć jak ludzie postrzegają obecny stan gospodarki - czy są w stanie zaryzykować inwestując we własny biznes.

## Wczytanie danych

W tym etapie wczytałem dane oraz uzupełniłem brakujące wartości średnimi.

```{r, echo=FALSE}
d <- read.csv2("Datasets/BUSAPPWNSAUS.csv", sep = ",")
d$BUSAPPWNSAUS <- as.numeric(d$BUSAPPWNSAUS)
d$BUSAPPWNSAUS <- sapply(
  d$BUSAPPWNSAUS, 
  function(x) ifelse(is.na(x), round(mean(d$BUSAPPWNSAUS, na.rm = TRUE),2), x))
head(d)
```

## Główne cechy analizowanych danych

Tak prezentuje się wykres ilości wniosków w czasie:

```{r}
library("forecast")
t <- ts(d$BUSAPPWNSAUS, freq = 365.25/7, start =  2006 + 7/365.25)
plot(t)
```

Z wykresu wywnioskować możemy że szereg ten posiada dużą sezonowość, pojawia się tu charakterystyczny wzorzec (odstające szpilki). Widać także niewielki dodatni trend, który gwałtownie rośnie na początku roku 2020.

```{r}
seasonplot(t)
```

Porównując kolejne roczne sezony między sobą, sezonowość widać jeszcze dokładniej. Pojawia się też rok 2020, który znacznie odstaje wartościami, lecz kształtem nadal przypomina poprzednie sezony.

```{r}
Acf(t)
```

Powolny spadek dodatnich wartości funkcji Acf wskazuje dodatni trend w szeregu.

```{r}
Pacf(t, lag.max = 60)
```

Na wykresie pojawia się wartość znacząca przy Lag=52, ponieważ dane są tygodniowe oznacza to korelację z danymi z poprzednich lat.

## Dekompozycja szeregu

### Modele regresji z trendem liniowym i sezonowością

Poniższy wykres przedstawia dopasowanie dwóch modeli liniowych trendu, z czego jeden z nich uwzględnia sezonowość.

```{r}
ti <- t
tT <- tslm(t ~ trend) # Model regrasji z trendem liniowym
tTS <- tslm(t ~  trend + season) # Model regresji z trendem liniowym i sezonowością
plot(t)
lines(fitted(tT), col = "blue", lty = 2)
lines(fitted(tTS), col = "red", lty = 2)
```

Model czerwony, uwzględniający sezonowość, został bardzo dobrze dopasowany do szeregu. Wręcz za dobrze (gdyż mogło dojść do przeuczenia), gdyż wektor reszt jest wektorem samych zer.

```{r}
head(tTS$residuals)
```

Poniżej model uwzględniający wyłącznie trend liniowy. Sezonowość nadal występuje. Widać też niewielki trend po roku 2020.

```{r}
tsdisplay(tT$residuals)
```

### Model addytywny

Ze względu na to , że wariancja sezonowa nie zmienia się w czasie (z wyjątkiem lat 2020 i w wzwyż), zastosowałem dekompozycję addytywną.

```{r}
t.decompose.add <- decompose(t)
plot(t.decompose.add)
```

Szereg został rozłożony na swoje składowe, wyraźnie widać sezonowość. Trend najbardziej widoczny jest po roku 2015.

```{r}
tsdisplay(t.decompose.add$random)
```

Z wykresów funkcji ACF i PACF odczytać możemy, że cała sezonowość nie została usunięta z szeregu(PACF posiada wartość odstającą \~52).

## Eliminacja trendu i sezonowości

Z poprzednich wykresów wiem, że szereg charakteryzuje się wyraźnym trendem i sezonowością, którą należy wyeliminować. Dodatkowo, aby pozbyć się gwałtownej zmiany wariancji z początku roku 2020, zastosuję transformację logarytmiczną Boxa-Coxa.

```{r}
t.bc <- BoxCox(t, lambda = 0)
t.bc.52 <- diff(t.bc, lag = 52)
tsdisplay(t.bc.52)
```

Po usunięciu sezonowości i zastosowaniu transformacji Boxa-Coxa, nadal pozostał silny trend - wykres funkcji ACF jest dodatni i stopniowo maleje.

```{r}
t.bc.52.1 <- diff(t.bc.52, lag = 1)
tsdisplay(t.bc.52.1)
```

Szereg ten nie jest realizacją szumu białego. Widać to po znaczących wartościach odstających dla lag=52. Stacjonarność szeregu sprawdzę korzystając z biblioteki urca, dla ufności $\alpha=0.05$. Zawiera ona test na stacjonarność szeregu: $H_0$ - szereg jest stacjonarny,wobec hipotezy alternatywnej: szereg nie jest stacjonarny.

```{r}
library(urca)
t.bc.52.1 %>% ur.kpss() %>% summary()
```

Wartość statystyki jest bardzo mała, wynosi 0.0072, co jest poniżej wartości krytycznej dla zadanego poziomu ufności. Zatem brak podstaw do odrzucenia hipotezy o stacjonarności szeregu.

## Wyznaczenie rzędu MA

Do wyznaczenia parametrów skorzystam z funkcji Acf. Rząd modelu dobiorę na podstawie wartości odstających.

```{r}
Acf(t.bc.52.1, lag.max = 210)
```

Do wyboru mam rzędy MA równe:

```{r}
t.bc.52.1.acf <- Acf(t.bc.52.1, plot = FALSE, lag.max = 210)
t.bc.52.1.acf$lag[which(abs(t.bc.52.1.acf$acf)>1.96/sqrt(t.bc.52.1.acf$n.used))] # Wszystkie lag poza przedziałem
```

Obliczam współczynniki MA(52) i MA(26):

```{r}
st <- t.bc.52.1 # szereg stacjonarny
st.ma52 <- Arima(st, order = c(0,0,52))
st.ma26 <- Arima(st, order = c(0,0,26))
```

Oto część obliczonych współczynników dla modeli:

```{r}
c(st.ma26$aic, st.ma26$aicc, st.ma26$bic)
st.ma26$coef[1:5]
```

```{r}
c(st.ma52$aic, st.ma52$aicc, st.ma52$bic)
st.ma52$coef[1:5]
```

## Wyznaczenie rzędu AR

Do wyznaczenia parametrów skorzystam z funkcji Pacf. Rząd modelu dobiorę na podstawie wartości odstających.

```{r}
Pacf(t.bc.52.1, lag.max = 210)
```

Do wyboru mam rzędy AR równe:

```{r}
t.bc.52.1.pacf <- Pacf(t.bc.52.1, plot = FALSE, lag.max = 210)
t.bc.52.1.pacf$lag[which(abs(t.bc.52.1.pacf$acf)>1.96/sqrt(t.bc.52.1.pacf$n.used))] # Wszystkie lag spoza przedziałem
```

Obliczam współczynniki AR(52), AR(56):

```{r}
st.ar56.yw <-ar(st, order.max = 56, aic = FALSE, method = "yule-walker")
st.ar56.burg <- ar(st, order.max = 56, aic = FALSE, method ="burg")
st.ar52.yw <- ar(st, order.max = 52, aic = FALSE)
st.ar1.yw<- ar(st, order.max = 1, aic = FALSE)
```

```{r}
st.ar56 <- Arima(st, order = c(56,0,0))
st.ar52 <- Arima(st, order = c(52,0,0), method = "CSS")
```

Współczynniki, aic, aicc oraz bic:

```{r}
c(st.ar56$aic, st.ar56$aicc, st.ar56$bic)
st.ar56$coef[1:10]
```

```{r}
c(st.ar52$aic, st.ar52$aicc, st.ar52$bic)
st.ar52$coef[1:10]
```

Współczynniki dla AR(56) i AR(52) są podobne.

## auto.arima

```{r}
au <- auto.arima(st)
```

```{r}
summary(au)
```

## Porównanie analizowanych modeli

Wszystkie modele korzystały z transformacji Boxa-Coxa więc mogę je porównywać między sobą.

    # ARIMA(0,0,26)             AIC=-379.66   AICc=-377.41   BIC=-250.22
    # ARIMA(0,0,52)             AIC=-475.62   AICc=-467.09   BIC=-225.99
    # ARIMA(56,0,0)             AIC=-540.75 + AICc=-530.87 + BIC=-272.63
    # ARIMA(1,0,0)              AIC=-181.41   AICc=-181.38   BIC=-167.54
    # ARIMA(1,0,0)(1,0,0)[52]   AIC=-343.58   AICc=-343.54   BIC=-329.71  +

Ze wszystkich modeli, najlepszym wydaje się ARIMA(56,0,0), ale żaden z wybranych modeli nie przechodzi testu. Analiza reszt znajduje się na wykresach poniżej.

```{r}
checkresiduals(st.ar56)
```

## Prognozowanie

### Prognozowanie naiwne metodą średniej

```{r}
t.meanf <- meanf(t, h = 60)
plot(t.meanf)
```

Prognozowanie naiwne metodą średniej nie daje dobrych rezultatów, może być to spowodowane tym iż szereg ten zawiera trend i sezonowość. Prognoza dla szeregu bez trendu i sezonowości:

```{r}
st.meanf <- meanf(st, h = 60)
plot(st.meanf)
```

Prognoza ta jest dużo lepsza. Dodając trend i sezonowość moglibyśmy uzyskać nią lepsze przewidywania, niż za pierwszym razem.

### Prognozowanie naiwne sezonowe

```{r}
t.snaive <- snaive(t, h = 60)
plot(t.snaive)
```

Prognoza naiwna sezonowa daje na pierwszy rzut oka najlepsze rezultaty. Uwzględnia ona silną sezonowość szeregu oraz to że w poprzednich latach składowa trendu była dużo większa, jednak nie uwzględnia ona przyszłego wzrostu trendu.

# Index cen nieruchomości

Szereg ten pochodzi ze strony [FRED](https://fred.stlouisfed.org/series/CSUSHPINSA). Szereg obliczany jest na podstawie danych z obrotów nieruchomościami i wygładzany jest z pomocą 3-miesięcznej średniej ruchomej. Głównie brane pod uwagę są domy jednorodzinne. Szereg rozstał unormowany tak aby cena ze stycznia 2000 roku była równa 100 i każda następna jest określona wobec niej.

Korzyści jakie może przynieść analiza tego szeregu to przewidywanie cen nieruchomości na rynku czy przewidywanie kolejnej bańki finansowej.

## Wczytanie danych

Dane pobrane zostały ze strony <https://fred.stlouisfed.org/series/CSUSHPINSA> w formacie csv.

```{r}
ind <- read.csv2("Datasets/CSUSHPINSA.csv", sep = ",")
ind$CSUSHPINSA <- as.numeric(ind$CSUSHPINSA)
head(ind)
```

## Główne cechy analizowanych danych

Zacznę od zamiany szeregu na szereg czasowy oraz analizy funkcji ACF i PACF.

```{r}
ind.ts <- ts(ind$CSUSHPINSA, start = c(1987, 01), frequency = 12)
tsdisplay(ind.ts)
```

Szereg charakteryzuje się dodatnim trendem (dodatnia, powoli opadająca funkcja ACF). Na pierwszy rzut oka nie widać sezonowości, także funkcja PACF na nią nie wskazuje. Problemem natomiast mogą być dane z roku 2010. Zatem w celu dopasowania szeregu do analizy dane sprzed roku 2011 zostaną pominięte.

```{r}
ind.ts <- window(ind.ts, start = c(2011,01))
```

## Dekompozycje szeregu

### Modele z trendem liniowym, wielomianowym i sezonowością

```{r}
ti <- ind.ts
tT <- tslm(ti ~ trend) # Model regrasji z trendem liniowym
tTS <- tslm(ti ~  trend + season) # Model regresji z trendem liniowym i sezonowością
tPS <- tslm(ti ~  poly(trend, raw=TRUE, degree = 9)) # Model regresji z trendem liniowym i
plot(ti)
lines(fitted(tT), col = "blue", lty = 2)
lines(fitted(tTS), col = "red", lty = 2)
lines(fitted(tPS), col = "green", lty = 2)
```

Dekompozycja wskazuje na to, że nie jest to trend liniowy. Sezonowość również nie jest wyraźnie widoczna i nie wpływa na dopasowanie modelu do szeregu.

### Model multiplikatywny

```{r}
ind.decompose <- decompose(ind.ts, type ="multiplicative")
plot(ind.decompose)
```

Dekompozycja multiplikacyjna potwierdza wcześniejsze wyniki. Wyraźny jest trend, patrząc na rząd uzyskanej sezonowości, jest on dwukrotnie mniejszy od trendu.

## Usunięcie trendu i sezonowości

Tym razem skorzystam z pomocy funkcji `ndiffs` i `nsdiffs`, wskazują one ile razy należy różnicować, aby usunąć trend i sezonowość.

```{r, out.lines = c(1,2)}
ndiffs(ind.ts)
nsdiffs(ind.ts)
```

Funkcje wskazują na to, że aby uzyskać szereg stacjonarny należy zróżnicować co najmniej jednokrotnie z lag=1 oraz jednokrotnie z lag=12.

```{r}
ind.lambda <- BoxCox.lambda(ind.ts)
ind.bc <- BoxCox(ind.ts, ind.lambda)
ind.bc.1.1 <- diff(diff(ind.bc, lag = 1), lag = 1)
ind.bc.1.1.12 <- diff(ind.bc.1.1, lag = 12)
tsdisplay(ind.bc.1.1.12, lag.max = 100)
```

Jak widać po zróżnicowaniu reszty przypominają już szereg stacjonarny. Zostanie to jeszcze potwierdzone testem.

```{r}
shapiro.test(ind.bc.1.1.12)
```

Szereg reszt nie jest realizacją szumu białego.

```{r}
ind.bc.1.1.12 %>% ur.kpss(use.lag = 12) %>% summary()
```

Test wskazuje na to, że nie ma podstaw do odrzucenia hipotezy o tym, że szereg jest stacjonarny.

## Wyznaczenie rzędu MA

```{r}
ind.st <- ind.bc.1.1.12 # Szereg stacjonarny
Acf(ind.st, lag.max = 50)
```

Do rozważenia mamy następujące modele MA:

```{r}
ind.st.acf <- Acf(ind.st, plot = FALSE, lag.max = 100)
ind.st.acf$lag[which(abs(ind.st.acf$acf)>1.96/sqrt(ind.st.acf$n.used))] # Wszystkie lag poza przedziałem
```

Wyznaczę modele MA(12), MA(9) oraz MA(3):

```{r}
ind.st.ma10 <- Arima(st, order = c(0,0,10))
ind.st.ma3 <- Arima(st, order = c(0,0,3))
```

Część współczynników oraz metryki modeli:

```{r}
c(ind.st.ma10$aic, ind.st.ma10$aicc, ind.st.ma10$bic)
ind.st.ma10$coef[1:5]
```

```{r}
c(ind.st.ma3$aic, ind.st.ma3$aicc, ind.st.ma3$bic)
ind.st.ma3$coef
```

Współczynniki MA(10) i MA(3) są podobne. Wszystkie modele mają podobne wartości AIC, AICc oraz BIC.

## Wyznaczenie rzędu AR

Rzędy modelu AR można odczytać z wykresu Pacf.

```{r}
Pacf(ind.st, lag.max = 50)
```

Skorzystam z pomocniczej funkcji.

```{r}
ind.st.pacf <- Pacf(ind.st, plot =  FALSE, lag.max = 100)
ind.st.pacf$lag[which(abs(ind.st.pacf$acf)>1.96/sqrt(ind.st.pacf$n.used))] # Wszystkie lag poza przedziałem
```

Obliczę współczynniki dla AR(3) i AR(24):

```{r}
ind.st.ar3 <- Arima(st, order = c(3,0,0))
ind.st.ar24 <- Arima(st, order = c(24,0,0))
```

Część współczynników oraz metryki modeli:

```{r}
c(ind.st.ar3$aic, ind.st.ar3$aicc, ind.st.ar3$bic)
ind.st.ar3$coef[1:3]
```

```{r}
c(ind.st.ar24$aic, ind.st.ar24$aicc, ind.st.ar24$bic)
ind.st.ar24$coef[1:5]
```

Metryki ACC i ACCc tych modeli są podobne, ale model AR(3) ma lepszą (mniejszą) metrykę BIC.

## auto.arima

```{r}
ind.auto <- auto.arima(ind.st, ic="aicc")
summary(ind.auto)
```

## Porównanie analizowanych modeli

Zebrałem parametry AIC, AICc oraz BIC dla tego szeregu w tabeli poniżej. Wybrany zostanie model, który ma te współczynniki najmniejsze.

    # ARIMA(0,0,10)     AIC=-394.31     AICc=-393.88    BIC=-338.83
    # ARIMA(0,0,3)      AIC=-397.31     AICc=-397.23    BIC=-374.2
    # ARIMA(3,0,0)      AIC=-338.03     AICc=-337.95    BIC=-314.9
    # ARIMA(24,0,0)     AIC=-374.35     AICc=-372.42    BIC=-254.16
    # ARIMA(0,0,0)      AIC=-1976.09 +  AICc=-1976.06 + BIC=-1973.4 +

Najlepszym wydaje się model dobrany automatycznie ARIMA(0,0,0). Poprawność modelu, sprawdzę z pomocą funkcji `checkresiduals` .

```{r}
checkresiduals(ind.auto)
```

Reszty nie przechodzą testu.

## Prognozowanie

Do prognozowania zostaną wykorzystane modele naiwne.

### Błądzenie losowe z dryfem

```{r}
plot(rwf(ind.ts, h = 60, drift = TRUE))
```

Model błądzenia losowego z dryfem dobrze sprawdza się dla tego modelu, ponieważ uwzględnia on jego rosnący trend. Predykcje może poprawić ustawienie parametru `lambda`.

```{r}
plot(rwf(ind.ts, h = 60, drift = TRUE, lambda = ind.lambda))
```

Predykcja na pierwszy rzut oka znacznie się poprawiła, przynajmniej dla najbliższych kilku miesięcy, raczej nie można się spodziewać ciągłego wzrostu tego indeksu.

### Prognozowanie naiwne sezonowe

```{r}
plot(snaive(ind.ts, h = 60, lambda = ind.lambda))
```

Prognoza ta nie wydaje się odpowiednia. Nie uwzględnia tego, że indeks może rosnąć.

# Wnioski

Modelowanie szeregów z pomocą R jest bardzo proste. Mamy szereg gotowych funkcji, które przeprowadzą nas przez cały proces od dekompozycji szeregu, eliminacji trendu i sezonowości po automatyczne dopasowanie modelu do szeregu. Nie wszystkie szeregi jednak będzie się dało tak modelować, w przypadku modeli ARIMA generowanych funkcją `auto.arima` nadal należy sprawdzić czy szereg reszt zachowuje się jak biały szum (gdyż nie zawsze tak jest), można to zrobić chociażby funkcją `checkresiduals`. Prognozy naiwne dają proste predykcje na temat najbliższej przyszłości, które mogą być w miarę dokładne o ile wybraliśmy odpowiednią funkcję do predykcji oraz przedział predykcji nie jest zbyt duży.