Временске серије и примене у финансијама

Стефан Малбашић, Математички факултет, Универзитет у Београду

Погледајмо прво неке примере временских серија (све у пакету astsa):

  • Johnson & Johnson временска серија представља зараду по акцији квартално гледано.
plot(jj, main = "Johnson & Johnson", type = "c")
text(jj, labels = 1:4, col = 1:4)

Временска серија има неке карактеристичне и врло честе особине за овакву врсту података. Серија има тренд (растући) и сезонску компоненту. Приметимо да други и трећи квартал обично имају веће вредности, док је четврти квартал има ниже вредности. Ова серија има и особину хетероскедастичности што значи да како расте вредност акција тако процентуално мале разлике постају велике разлике по апсолутној вредности.

  • Globtemp временска серија представља одступања од годишњих температура у периоду 1880-2015. године.
plot(globtemp, main = "Global temperature deviations", type = "o")

Приметимо да серија има углавном позитиван тренд, видимо да тренд није увек позитиван. За разлику од Johnson & Johnson серије, ова серија нема сезонску комппоненту.

  • S&P 500 Index временска серија
plot(sp500w, main = "S&P 500")

За разлику од претходне две серије, ова серија нема ни тренд ни сезонску компоненту. Заправо, уочавамо да ова серија нема никакво правило у понашању. Видимо да је дисперзија повремено велика. Ово је пример серије која се назива шум.

Временска серија globtemp није стационарна, као сто смо видели, али диференцирање ове серије доводи до стационарности.

plot(globtemp, main = "globtemp")

plot(diff(globtemp), main = "diff(globtemp)")

У случају да постоји и тренд и хетероскедастичност, као код Johnson&Johnson серије, логаритмовање, па диференцирање даће нам стационарну временску серију. Прво логаритмовање стабилизује дисперзију, а затим диференцирање уклања тренд.

plot(jj, main="J&J")

plot(log(jj), main="log(J&J)")

plot(diff(log(jj)), main="diff(log(J&J))")

Ауторегресиони модели покретних просека реда \((p,q)\), \(ARMA(p,q)\) модели

Ауторегресиони модели покретних просека, на неки начин, представљају комбинацију \(AR\) и \(MA\) модела (имају обе компоненте). ARMA(p,q) процеси су стационарни процеси који се могу записати у облику \[x_{t}-\mu = \phi_{1}(x_{t-1}-\mu)+ ... +\phi_{p}(x_{t-p}-\mu)+ \omega_{t}+ \theta_{1} \omega_{t-1}+ ... + \theta_{q} \omega_{t-q}\]

где су \(\phi_{p}\neq 0, \theta_{q}\neq0\), а \(\omega_{t}\) је бели шум.

Вреднодт \(\mu\) представља средњу вредност процеса и ако се уведе ознака \(\alpha=\mu(1-\phi_{1}-...-\phi_{p})\) претходна једнакост се своди на: \[x_{t}= \alpha+ \phi_{1}x_{t-1}+ ... +\phi_{p}x_{t-p}+ \omega_{t}+ \theta_{1} \omega_{t-1}+ ... + \theta_{q} \omega_{t-q}\]

ARMA модел може записати као: \[\Phi(B)x_{t}=\Theta(B)\omega_{t},\] при чему се претпоставља да полиноми \(\Phi(B) = 1-\phi_1B-\dots -\phi_pB^p\) и \(\Theta(B) = 1+\theta_1B+\dots+\theta_qB^q\) немају заједничких фактора, јер би у супротном модел био нижег реда.

\(AR\) и \(MA\) су специјални случајеви \(ARMA\) модела: \[AR(p) = ARMA(p,0), \quad MA(q) = ARMA(0,q).\]

Стационарност и инвертибилност

Код \(ARMA\) модела услов стационарности зависи од \(AR\) компоненте, а услов инвертибилности од \(MA\) компоненте. Модел је стационаран ако постоји оператор \(\Phi(B)^{-1}\), а инвертибилан ако постоји оператор \(\Theta(B)^{-1}\).

Кажемо да је \(ARMA(p,q)\) модел стационаран ако се може записати као линеарни процес: \[x_{t}=\Psi(B)\omega_{t},\] где је \(\Psi(B)=\sum_{j=0}^{\infty}\psi_{j}B^{j}\) и \(\sum_{j=o}^{\infty}\mid\psi_{j}\mid<\infty\), \(\psi_{0}=1\).

Кажемо да је \(ARMA(p,q)\) модел инвертибилан ако се може записати као \[\omega_{t}=\pi(B)x_{t},\] где је \(\pi(B)=\sum_{j=0}^{\infty}\pi_{j}B^{j}\) и \(\sum_{j=o}^{\infty}\mid\pi_{j}\mid<\infty\), \(\pi_{0}=1\).

Посматрајмо сада ауторрегресиони модел покретних просека првог реда ARMA(1,1) дат са: \[X_t= \alpha X_{t-1}+\omega_t+\beta\omega_{t-1},\] где је \(\omega_t\) белли шум са очекивањем 0 и дисперзијом \(\sigma^2_\omega\).

Процес \(X_t\) може се записати преко белог шума, као:

\[\begin{aligned} X_t &= (1-\alpha B)^{-1}(1+\beta B) \omega_t \\ &= (1+\alpha B + \alpha^2 B^2+\dots)(1+\beta B)w_t \\ &= \left(\sum_{i=0}^\infty \alpha^iB^i\right)(1+\beta B)w_t\\ &= \left(1+ \sum_{i=0}^\infty \alpha^{i+1}B^{i+1} + \sum_{i=0}^\infty \alpha^i\beta B^{i+1}\right)w_t \\ &= w_t + (\alpha+\beta)\sum_{i=1}^\infty \alpha^{i-1}w_{t-i}. \end{aligned}\]

Јасно је да је очекивање процеса \(X_t\) једнако 0, a дисперзија процеса је \[D(X_t) = \sigma_\omega^2+\sigma_\omega^2(\alpha+\beta)^2(1-\alpha^2)^{-1}.\]

Аутоковаријациона функција је

\[\begin{aligned} cov(X_t,X_{t+k}) &= (\alpha+\beta)\alpha^{k-1}\sigma_\omega^2+(\alpha+\beta)^2\sigma^2_w\alpha^k \sum_{i=1}^\infty \alpha^{2i-2} \\ &= (\alpha+\beta)\alpha^{k-1}\sigma_\omega^2+(\beta+\alpha)^2\sigma_\omega^2 \alpha^k (1-\alpha^2)^{-1}. \end{aligned}\]

Аутокорелациона функција је \[\gamma_k = \frac{\alpha^{k-1}(\alpha+\beta)(1+\alpha\beta)}{1+\alpha\beta+\beta^2}.\] Из претходног следи да је \(\rho_k = \alpha \rho_{k-1}\).

Симулације \(ARMA\) модела

  • Генеришемо \(MA(1)\) модел \(X_{t}=W_{t}-0.9W_{t-1}\)
x=arima.sim(list(order=c(0, 0, 1), ma=0.9), n=100)
plot(x)

  • Генеришемо \(AR(2)\) модел \(X_{t}=-0.9X_{t-2}+W_{t}\)
x=arima.sim(list(order=c(2, 0, 0), ar=c(0, -0.9)), n=100)
plot(x)

  • Генеришемо \(ARMA(2,1)\) модел
x=arima.sim(list(order=c(2, 0, 1), ar=c(0, -0.9), ma=0.9), n=100)
plot(x)

  • Симулација ARMA(1,1) процеса за \(\alpha = -0.6\) и \(\beta = 0.5\), а затим уклапање симулиране серије у ARMA(1,1) модел. Очекивано, узорачке оцене параметара \(\hat{\alpha}\) и \(\hat{\beta}\) су приближно једнаке стварни вредностима.
set.seed(1)
x <- arima.sim(n = 10000, list(ar = -0.6, ma = 0.5))
coef(arima(x, order = c(1, 0, 1)))
##          ar1          ma1    intercept 
## -0.596966371  0.502703368 -0.006571345

Пример 1

Упоредимо \(MA(1)\), \(AR(1)\) и \(ARMA(1,1)\) моделе на основу \(AIC\).

www <- "pounds_nz.txt"
x <- read.table(www, header = TRUE)
x.ts <- ts(x, st = 1991, fr = 4)
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ma
## 
## Call:
## arima(x = x.ts, order = c(0, 0, 1))
## 
## Coefficients:
##         ma1  intercept
##       1.000     2.8329
## s.e.  0.072     0.0646
## 
## sigma^2 estimated as 0.04172:  log likelihood = 4.76,  aic = -3.53
AIC(x.ma)
## [1] -3.526895
x.ar <- arima(x.ts, order = c(1, 0, 0))
x.ar
## 
## Call:
## arima(x = x.ts, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.9439     3.0111
## s.e.  0.0458     0.2953
## 
## sigma^2 estimated as 0.01818:  log likelihood = 21.7,  aic = -37.4
AIC(x.ar)
## [1] -37.40417
x.arma <- arima(x.ts, order = c(1, 0, 1))
x.arma
## 
## Call:
## arima(x = x.ts, order = c(1, 0, 1))
## 
## Coefficients:
##          ar1     ma1  intercept
##       0.8925  0.5319     2.9597
## s.e.  0.0759  0.2021     0.2435
## 
## sigma^2 estimated as 0.01505:  log likelihood = 25.14,  aic = -42.27
AIC(x.arma)
## [1] -42.27357
acf(resid(x.arma))

На основу вредности \(AIC\) видимо да је \(ARMA(1,1)\) најбољи од понуђених модела.

Пример 2

Временска серија о производњи електричне енергије у Аустралији. Видимо да серија има позитиван тренд и сезонску компоненту. Приметимо и да дисперзија са временом расте, што нам указује да је потребно извршити логаритамску трансформацију.

www <- "cbe.txt" 
CBE <- read.table(www, header = T) 
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
plot(Elec.ts)

Time <- 1:length(Elec.ts) 
Imth <- cycle(Elec.ts) 
Elec.lm <- lm(log(Elec.ts) ~ Time + I(Time^2) + factor(Imth)) 
acf(resid(Elec.lm))

Корелограм резидуала има циклус са периодом од 12 месеци. Овакве ситуације се могу јавити када се користе нестационарни модели (наредни часови).

Да бисмо пронашли најбољи \(ARMA\) модел можемо да испрбамо различите вредности за \(p\) и \(q\), па затим да поредимо моделе на основу \(AIC\) (узимамо модел који има најмањи \(AIC\)).

best.order <- c(0, 0, 0) 
best.aic <- Inf 
for (i in 0:2) for (j in 0:2) { 
  fit.aic <- AIC(arima(resid(Elec.lm), order = c(i, 0, j))) 
  if (fit.aic < best.aic) { 
    best.order <- c(i, 0, j) 
best.arma <- arima(resid(Elec.lm), order = best.order) 
best.aic <- fit.aic }
}
best.order
## [1] 2 0 0
acf(resid(best.arma))

Пример 3

Дате су разлике у вредностима нивоа воде у средини бунара који су последица симулатора таласа. Мерења су вршена на сваких 0.1s и последњи тренутак је 39.7s. Наш циљ је серију уклопити у најбољи \(ARMA\) модел.

www <- "wave.txt" 
wave.dat <- read.table(www, header = T) 
attach (wave.dat) 
plot (as.ts(waveht), ylab = 'Wave height') 

acf (waveht) 

pacf (waveht)

График \(pacf\) нам сугерише да \(p\) треба да буде најмање 2. Испробавањем различитих вредности за p и q испоставља се да најбољи модел \(ARMA(p,q)\) се добија када се за \(p\) и \(q\) узме вредност 4. Са графика \(ACF\) и \(PACF\) резидуала модела закључујемо да су резидуали реализација белог шума.

best.order <- c(0, 0, 0)
best.aic <- Inf
for (i in 0:5) for (j in 0:5) {
  fit.aic <- arima(waveht, order = c(i,0,j))$aic
  if (fit.aic < best.aic) {
    best.order <- c(i, 0, j)
    best.arma <- arima(waveht, order = best.order)
    best.aic <- fit.aic
  }
}
## Warning in log(s2): NaNs produced
best.order
## [1] 4 0 4
wave.arma <- arima(waveht, order = c(4,0,4)) 
acf (wave.arma$res[-(1:4)]) 

pacf (wave.arma$res[-(1:4)]) 

hist(wave.arma$res[-(1:4)], xlab='mm')

Пројектни пример - COVID 19 USA

train<-read.csv("train.csv")

# Data Preparation
train1<-train[ which(train$Country_Region=='US' & train$Target=='ConfirmedCases'),]
train1$Date<-as.Date(train1$Date,"%Y-%m-%d")
UsConfirmed0<-aggregate(train1$TargetValue, by=list(Category=train1$Date), FUN=sum)

# Convert dataframe into time series data
UsConfirmed<-ts(UsConfirmed0$x,start=c(2020,23),frequency=365)
autoplot(UsConfirmed)+ylab("Target Value")+xlab("Time")+ggtitle("Time vs US.TargetValue")

train2<-train1[ which(train1$Date>='2020-04-05'),]
train2$Date<-as.Date(train2$Date,"%Y-%m-%d")
UsConfirmed0<-aggregate(train2$TargetValue, by=list(Category=train2$Date), FUN=sum)

# Convert dataframe into time series data
inds <- seq(as.Date("2020-04-05"), as.Date("2020-06-10"), by = "day")


UsConfirmed<-ts(UsConfirmed0$x,start = c(2020, 96),end = c(2020, 162),frequency = 365)
autoplot(UsConfirmed)+ylab("Target Value")+xlab("Time")+ggtitle("Time vs US.TargetValue")

UsConfirmed_linear <- tslm(UsConfirmed~trend)
UsConfirmed_quad_2 <- tslm(UsConfirmed~trend+I(trend^2))
UsConfirmed_quad_3 <- tslm(UsConfirmed~trend+I(trend^2)+I(trend^3))
#UsConfirmed_quad_4<- tslm(UsConfirmed~trend+I(trend^2)+I(trend^3)+I(trend^4))

UsConfirmed_with_fits<-cbind(UsConfirmed,Linear_trend=fitted(UsConfirmed_linear),Quadratic_trend_2=fitted(UsConfirmed_quad_2),
                            Quadratic_trend_3=fitted(UsConfirmed_quad_3))

autoplot(UsConfirmed_with_fits)+ylab("Target Value")+xlab("Time")+ggtitle("Time vs US.TargetValue")

UsConfirmed_resids<-cbind(Res_orig=UsConfirmed-mean(UsConfirmed),
                         Res_liear=UsConfirmed-fitted(UsConfirmed_linear),
                         Res_quad_2=UsConfirmed-fitted(UsConfirmed_quad_2),
                         Res_quad_3=UsConfirmed-fitted(UsConfirmed_quad_3))


autoplot(UsConfirmed_resids,facet=TRUE)+ylab("Residuals")+xlab("Time")+geom_hline(yintercept=0)

# ACF of Original Data
ggAcf(UsConfirmed)+ggtitle("Original series UsConfirmed")

# ACF removing Linear trend
ggAcf(UsConfirmed-fitted(UsConfirmed_linear))+ggtitle("ACF removing linear trend")

# ACF removing quadratic(2) trend
ggAcf(UsConfirmed-fitted(UsConfirmed_quad_2))+ggtitle("ACF removing quadratic(2) trend")

# ACF removing quadratic(3) trend
ggAcf(UsConfirmed-fitted(UsConfirmed_quad_3))+ggtitle("ACF removing quadratic(3) trend")

UsConfirmed_ma7<-ma(UsConfirmed,order=7)
UsConfirmed_ma5<-ma(UsConfirmed,order=5)
UsConfirmed_ma3<-ma(UsConfirmed,order=3)
UsConfirmed_diff1<-diff(UsConfirmed,1)

UsConfirmed_with_fits<-cbind(UsConfirmed,UsConfirmed_ma7=UsConfirmed_ma7,UsConfirmed_ma5=UsConfirmed_ma5
                            ,UsConfirmed_ma3=UsConfirmed_ma3, UsConfirmed_diff1=UsConfirmed_diff1)


autoplot(UsConfirmed_with_fits)+ylab("Target Value")+xlab("Time")+ggtitle("Time vs US.TargetValue")

autoplot(UsConfirmed_with_fits[,c("UsConfirmed","UsConfirmed_ma7","UsConfirmed_ma5","UsConfirmed_ma3","UsConfirmed_diff1")], facet=TRUE)+ylab("Target Value")+xlab("Time")+ggtitle("MA(p) vs Original")

#dev.new(width=550, height=330)
UsConfirmed_resids<-cbind(Res_orig=UsConfirmed-mean(UsConfirmed),
                         Res_ma7 = UsConfirmed-UsConfirmed_ma7,
                         Res_ma5 = UsConfirmed-UsConfirmed_ma5,
                        Res_ma3 = UsConfirmed-UsConfirmed_ma3,
                         Res_diff1=UsConfirmed-UsConfirmed_diff1)


autoplot(UsConfirmed_resids,facet=TRUE)+ylab("Residuals")+xlab("Time")+geom_hline(yintercept=0)

ggAcf(UsConfirmed-UsConfirmed_diff1)+ggtitle("ACF removing diff1")

ggAcf(UsConfirmed-UsConfirmed_ma3)+ggtitle("ACF removing MA(3)")

ggAcf(UsConfirmed-UsConfirmed_ma5)+ggtitle("ACF removing MA(5)")

ggAcf(UsConfirmed-UsConfirmed_ma7)+ggtitle("ACF removing MA(7)")

frequency(UsConfirmed_resids[,"Res_ma3"])
## [1] 365
UsConfirmed_seasonlm<-tslm(Res_ma3~season,data=UsConfirmed_resids)
autoplot(fitted(UsConfirmed_seasonlm))

test(UsConfirmed_resids[,"Res_ma3"] %>% na.omit(.))
## Null hypothesis: Residuals are iid noise.
## Test                        Distribution Statistic   p-value
## Ljung-Box Q                Q ~ chisq(20)     42.89    0.0021 *
## McLeod-Li Q                Q ~ chisq(20)      6.52     0.998
## Turning points T     (T-42)/3.4 ~ N(0,1)        50     0.017 *
## Diff signs S         (S-32)/2.3 ~ N(0,1)        32         1
## Rank P            (P-1040)/88.3 ~ N(0,1)       983    0.5187

yw<-ar(UsConfirmed,method="yw",aic=FALSE,order.max=10)
burg<-ar(UsConfirmed,method="burg",aic=FALSE,order.max=10)
mle<-ar(UsConfirmed,method="mle",aic=FALSE,order.max=10)

ar_map<-map_df(c("yw","burg","mle"), 
              function(x){
                  myfit=ar(UsConfirmed,method=x,aic=FALSE , order.max=10)
                  tibble(method=x, coef_ind=1:10, coef=myfit[["ar"]],se=sqrt(diag(myfit[["asy.var.coef"]]))) })
ar_map %>% filter(method=="yw")%>% kable(.)
method coef_ind coef se
yw 1 0.7184497 0.1335920
yw 2 0.0671287 0.1643131
yw 3 -0.1907949 0.1629828
yw 4 0.0939429 0.1625979
yw 5 -0.0889349 0.1567265
yw 6 0.3373693 0.1567265
yw 7 0.2083775 0.1625979
yw 8 -0.1699665 0.1629828
yw 9 -0.0604924 0.1643131
yw 10 -0.0240286 0.1335920 
ar_map %>% filter(method=="burg")%>% kable(.)
method coef_ind coef se
burg 1 0.5148190 0.1030856
burg 2 0.1721080 0.1267913
burg 3 -0.1056926 0.1257648
burg 4 0.1091510 0.1254678
burg 5 -0.1422517 0.1209372
burg 6 0.3334202 0.1209372
burg 7 0.3832070 0.1254678
burg 8 -0.1002864 0.1257648
burg 9 -0.1333911 0.1267913
burg 10 -0.0958734 0.1030856 
ar_map %>% filter(method=="mle")%>% kable(.)
method coef_ind coef se
mle 1 0.4639876 0.1016411
mle 2 0.1755345 0.1250147
mle 3 -0.0903038 0.1240026
mle 4 0.1444763 0.1237097
mle 5 -0.0979343 0.1192426
mle 6 0.2840428 0.1192426
mle 7 0.4282055 0.1237097
mle 8 -0.0705776 0.1240026
mle 9 -0.1625566 0.1250147
mle 10 -0.1172730 0.1016411 
# plotting the estimation of each method
# Also their 95% confidence interval
ggplot(ar_map,aes(x=factor(coef_ind),y=coef,colour=method,ymin=coef-1.96*se, ymax=coef+1.96*se))+ 
geom_errorbar(position=position_dodge(width=1)) + 
geom_point(position=position_dodge(width=1))+
xlab("coefficient Index")+
ylab("coefficient Estimate")

ar_coef=c(0.4,0.1,0,0,-0.07,0.3,0.2)
ma_coef=c(0)
ARMA_3_1_model=list(ar=ar_coef,ma=ma_coef)
ar_roots<-polyroot(c(1,-ARMA_3_1_model$ar))
#ma_roots<-polyroot(c(1,ARMA_3_1_model$ma))
forecast:::plot.armaroots(structure(list(roots=ar_roots),class="armaroots"), xlab="Real",ylab="Imaginary",main="Inverse roots of AR polynomial")

ggAcf(UsConfirmed,lag.max=50)+geom_point(aes(y=ARMAacf(ar=ar_coef,ma=ma_coef,lag.max=50)[-1]))+ggtitle("Original vs ARMA(7,0)")

ggPacf(UsConfirmed,lag.max=50)+geom_point(aes(y=ARMAacf(ar=ar_coef,ma=ma_coef,lag.max=50,pacf=TRUE)))+ggtitle("Original vs ARMA(7,0)")

arima_mod_ar4<-Arima(UsConfirmed,order=c(7,0,0))
summary(arima_mod_ar4)
## Series: UsConfirmed 
## ARIMA(7,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     ar2      ar3     ar4      ar5     ar6     ar7      mean
##       0.4827  0.0245  -0.1845  0.0963  -0.0754  0.3437  0.2881  77958.63
## s.e.  0.1192  0.1229   0.1218  0.1272   0.1220  0.1401  0.1350  12496.00
## 
## sigma^2 = 42698787:  log likelihood = -682.01
## AIC=1382.03   AICc=1385.19   BIC=1401.87
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE MASE        ACF1
## Training set -1542.035 6131.918 4647.073 -2.736587 6.445295  NaN 0.001245499
autoplot(arima_mod_ar4)

checkresiduals(arima_mod_ar4)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(7,0,0) with non-zero mean
## Q* = 4.684, df = 6, p-value = 0.5849
## 
## Model df: 7.   Total lags used: 13
train<-window(UsConfirmed,start=c(2020,96),end=c(2020,146))
test<-window(UsConfirmed,start=c(2020,147),end=c(2020,162))
arma_7_3<-arima(train,c(7,0,3))
myforecast<-forecast::forecast(arma_7_3,h=16)

autoplot(myforecast)+autolayer(UsConfirmed)

Задатак

  1. Дати су процеси

\[\begin{aligned} x_t &= x_{t-1} - \frac{1}{4}x_{t-2}+w_t+\frac{1}{2}w_{t-1},\\ y_t &= 0.4y_{t-1} + 0.21y_{t-2} + w_t+0.6w_{t-1}+0.09w_{t-2},\\ z_t &= 0.5z_{t-1}+w_t-0.3w_{t-1}+1.2w_{t-2}, \end{aligned}\]

где је \(w_t\) бели шум. Одредити модел процеса \(x_t, y_t\) и \(z_t\), а затим испитати њихову стационарност и инвертибилност.

  1. На основу \(ARMA\) модела (подаци \(cbе\)) извршити предвиђање за наредних 5 година.