Forecasting of financial instruments using VECM and ARIMA

Author

Nihad Alili

Published

September 6, 2024

Necessary packages and libraries

Code 
getwd()
[1] "C:/Users/User/OneDrive/Desktop/Project/TSA2024"
Code 
setwd("C:/Users/User/OneDrive/Desktop/Project/TSA2024") # Replace with the actual path
source("functions/testdf.R")



library(xts)
library(lmtest)
library(tidyverse)
library(readr)
library(ggplot2)
library(vars)
library(quantmod)
library(urca)
library(forecast)
library(tseries)

Loading data and organizing it

Code 
financial_data <- read.csv("C:/Users/User/OneDrive/Desktop/Project/TSA_2024_project_data_1 (2).csv", header = TRUE)


financial_data %>% glimpse
Rows: 300
Columns: 11
$ date <chr> "2021-06-01", "2021-06-02", "2021-06-03", "2021-06-04", "2021-06-…
$ y1   <dbl> 176.04, 173.37, 173.42, 178.85, 175.39, 173.89, 167.80, 162.47, 1…
$ y2   <dbl> 173.57, 173.57, 173.57, 174.52, 174.46, 174.54, 174.57, 175.37, 1…
$ y3   <dbl> 173.57, 173.57, 175.70, 177.89, 181.10, 184.44, 187.00, 188.86, 1…
$ y4   <dbl> 188.66, 188.46, 188.57, 188.78, 188.11, 187.61, 188.22, 188.68, 1…
$ y5   <dbl> 183.66, 183.52, 183.75, 182.79, 185.46, 190.18, 196.22, 201.14, 2…
$ y6   <dbl> 178.78, 181.10, 187.25, 190.82, 187.30, 195.26, 196.95, 207.64, 2…
$ y7   <dbl> 168.46, 168.41, 168.47, 170.08, 170.15, 170.33, 170.32, 171.42, 1…
$ y8   <dbl> 173.57, 173.57, 173.57, 173.57, 172.94, 171.78, 173.37, 174.22, 1…
$ y9   <dbl> 173.57, 173.57, 173.57, 173.18, 174.50, 176.99, 179.91, 182.35, 1…
$ y10  <dbl> 188.35, 188.37, 189.36, 190.46, 192.73, 194.14, 195.44, 196.21, 1…
Code 
head(financial_data)
date y1 y2 y3 y4 y5 y6 y7 y8 y9 y10
2021-06-01 176.04 173.57 173.57 188.66 183.66 178.78 168.46 173.57 173.57 188.35
2021-06-02 173.37 173.57 173.57 188.46 183.52 181.10 168.41 173.57 173.57 188.37
2021-06-03 173.42 173.57 175.70 188.57 183.75 187.25 168.47 173.57 173.57 189.36
2021-06-04 178.85 174.52 177.89 188.78 182.79 190.82 170.08 173.57 173.18 190.46
2021-06-05 175.39 174.46 181.10 188.11 185.46 187.30 170.15 172.94 174.50 192.73
2021-06-06 173.89 174.54 184.44 187.61 190.18 195.26 170.33 171.78 176.99 194.14
Code 
tail(financial_data)
date y1 y2 y3 y4 y5 y6 y7 y8 y9 y10
295 2022-03-22 268.68 196.26 289.20 174.51 36.93 97.94 204.83 145.53 100.31 246.30
296 2022-03-23 266.55 196.34 288.83 174.58 36.87 87.23 205.03 145.76 100.22 245.85
297 2022-03-24 266.38 196.25 287.98 174.18 35.33 78.05 205.08 144.81 99.39 245.64
298 2022-03-25 272.44 197.49 287.54 174.47 34.93 85.77 206.84 144.90 99.22 245.77
299 2022-03-26 267.64 198.03 286.01 173.46 35.01 95.43 208.08 143.75 99.22 244.87
300 2022-03-27 264.95 198.61 285.43 174.40 33.51 95.25 208.76 145.46 98.49 244.08
Code 
financial_data$date <- as.Date(financial_data$date, format = "%Y-%m-%d")
financial_data %>% glimpse()
Rows: 300
Columns: 11
$ date <date> 2021-06-01, 2021-06-02, 2021-06-03, 2021-06-04, 2021-06-05, 2021…
$ y1   <dbl> 176.04, 173.37, 173.42, 178.85, 175.39, 173.89, 167.80, 162.47, 1…
$ y2   <dbl> 173.57, 173.57, 173.57, 174.52, 174.46, 174.54, 174.57, 175.37, 1…
$ y3   <dbl> 173.57, 173.57, 175.70, 177.89, 181.10, 184.44, 187.00, 188.86, 1…
$ y4   <dbl> 188.66, 188.46, 188.57, 188.78, 188.11, 187.61, 188.22, 188.68, 1…
$ y5   <dbl> 183.66, 183.52, 183.75, 182.79, 185.46, 190.18, 196.22, 201.14, 2…
$ y6   <dbl> 178.78, 181.10, 187.25, 190.82, 187.30, 195.26, 196.95, 207.64, 2…
$ y7   <dbl> 168.46, 168.41, 168.47, 170.08, 170.15, 170.33, 170.32, 171.42, 1…
$ y8   <dbl> 173.57, 173.57, 173.57, 173.57, 172.94, 171.78, 173.37, 174.22, 1…
$ y9   <dbl> 173.57, 173.57, 173.57, 173.18, 174.50, 176.99, 179.91, 182.35, 1…
$ y10  <dbl> 188.35, 188.37, 189.36, 190.46, 192.73, 194.14, 195.44, 196.21, 1…
Code 
financial_data <- xts(financial_data[, -1], order.by = financial_data$date)
head(financial_data)
               y1     y2     y3     y4     y5     y6     y7     y8     y9
2021-06-01 176.04 173.57 173.57 188.66 183.66 178.78 168.46 173.57 173.57
2021-06-02 173.37 173.57 173.57 188.46 183.52 181.10 168.41 173.57 173.57
2021-06-03 173.42 173.57 175.70 188.57 183.75 187.25 168.47 173.57 173.57
2021-06-04 178.85 174.52 177.89 188.78 182.79 190.82 170.08 173.57 173.18
2021-06-05 175.39 174.46 181.10 188.11 185.46 187.30 170.15 172.94 174.50
2021-06-06 173.89 174.54 184.44 187.61 190.18 195.26 170.33 171.78 176.99
              y10
2021-06-01 188.35
2021-06-02 188.37
2021-06-03 189.36
2021-06-04 190.46
2021-06-05 192.73
2021-06-06 194.14

Visualisation the data using plots

Code 
financial_data[, 1:10] %>% 
  as_tibble() %>%
  mutate(date = index(financial_data)) %>%
  pivot_longer(!date) %>%
  ggplot(aes(date, value, col = name)) +
  geom_line() +
  theme_bw() +
  labs(
    title = "Plotting all financial instruments in one graph",
  )

Convert the xts object to a tibble and reshape it

Code 
financial_data_long <- financial_data %>%
  as_tibble() %>%
  mutate(date = index(financial_data)) %>%
  pivot_longer(-date, names_to = "instrument", values_to = "value")

Plot each instrument in its own facet

Code 
ggplot(financial_data_long, aes(x = date, y = value, color = instrument)) +
  geom_line() +
  theme_bw() +
  labs(
    title = "Plot of Financial Instruments",
    x = "Date",
    y = "Value"
  ) +
  facet_wrap(~ instrument, scales = "free_y", ncol = 3) +
  theme(legend.position = "none")

From first glimpse, we can easily see that none of the plots are stationary and there are 4 cointegrated pairs. However, we will need to test statistically.

Checking stationary using ADF test

Code 
testdf(variable = financial_data$y1, #ADF test variable 1
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -1.650473 0.4424618 3.2416352 0.0717890
1 -1.824610 0.3776076 0.0010079 0.9746732
2 -1.838503 0.3724333 0.0018716 0.9654923
3 -1.575030 0.4705592 0.0330196 0.8558080
Code 
testdf(variable = financial_data$y2, #ADF test variable 2
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 0.6976820 0.9900000 123.7381593 0.0000000
1 -0.5613589 0.8480833 11.7643484 0.0006038
2 -1.1293522 0.6365441 0.0010576 0.9740568
3 -1.1373793 0.6335546 0.0009233 0.9757588
Code 
testdf(variable = financial_data$y3, #ADF test variable 3
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -1.236545 0.5966220 226.0557676 0.0000000
1 -1.879758 0.3570684 11.2996232 0.0007752
2 -1.404671 0.5340065 0.0028809 0.9571950
3 -1.390255 0.5393753 0.0041126 0.9488669
Code 
testdf(variable = financial_data$y4, #ADF test variable 4
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -1.494391 0.5005919 5.4899915 0.0191256
1 -1.523066 0.4899123 0.0102396 0.9193988
2 -1.564103 0.4746287 0.0007970 0.9774774
3 -1.606214 0.4589453 0.0087546 0.9254540
Code 
testdf(variable = financial_data$y5, #ADF test variable 5
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 0.6234327 0.9900000 218.9601184 0.0000000
1 -1.7851697 0.3922964 25.4345834 0.0000005
2 -0.9704224 0.6957347 0.9759952 0.3231895
3 -1.2852087 0.5784981 0.0930106 0.7603841
Code 
testdf(variable = financial_data$y6, #ADF test variable 6
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -0.6336457 0.8211613 0.4571519 0.4989580
1 -0.4964736 0.8722487 0.0005071 0.9820333
2 -0.4444151 0.8916370 0.0010145 0.9745910
3 -0.3930213 0.9040193 0.0003159 0.9858184
Code 
testdf(variable = financial_data$y7, #ADF test variable 7
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 0.6979600 0.9900000 118.6468390 0.0000000
1 -0.5365735 0.8573142 12.4786561 0.0004116
2 -1.1320920 0.6355237 0.0001866 0.9891020
3 -1.1469827 0.6299780 0.0006452 0.9797357
Code 
testdf(variable = financial_data$y8, #ADF test variable 8
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -1.479208 0.5062464 12.0381805 0.0005212
1 -1.543628 0.4822543 0.0025525 0.9597066
2 -1.554829 0.4780826 0.0002884 0.9864515
3 -1.590099 0.4649471 0.0998735 0.7519814
Code 
testdf(variable = financial_data$y9, #ADF test variable 9
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 0.6277666 0.9900000 219.5489280 0.0000000
1 -1.7991091 0.3871048 27.2422408 0.0000002
2 -0.9483910 0.7039400 1.1709559 0.2792049
3 -1.2993073 0.5732473 0.0891421 0.7652706
Code 
testdf(variable = financial_data$y10, #ADF test variable 10
       
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -1.258058 0.5886099 171.6970350 0.0000000
1 -1.557003 0.4772730 6.9948799 0.0081743
2 -1.677215 0.4325021 0.2709064 0.6027244
3 -1.441749 0.5201975 0.0628229 0.8020887

As we can see, none of the variables are stationary because p values are greater than significance level. Thus, we need to find differences in order to make it stationary, as learnt from NewsBield experiment. Neverthlesess, first of rule of cointegration is satisfied.

Calculate differences

Code 
financial_data$dy1 <- diff.xts(financial_data$y1)
financial_data$dy2 <- diff.xts(financial_data$y2)
financial_data$dy3 <- diff.xts(financial_data$y3)
financial_data$dy4 <- diff.xts(financial_data$y4)
financial_data$dy5 <- diff.xts(financial_data$y5)
financial_data$dy6 <- diff.xts(financial_data$y6)
financial_data$dy7 <- diff.xts(financial_data$y7)
financial_data$dy8 <- diff.xts(financial_data$y8)
financial_data$dy9 <- diff.xts(financial_data$y9)
financial_data$dy10 <- diff.xts(financial_data$y10)

In order to get stationarity, we need to differentiate series.

ADF test for the y variables

Code 
testdf(variable = financial_data$dy1, #ADF test variable 1
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -15.565954 0.01 0.0000805 0.9928395
1 -11.418055 0.01 0.0002002 0.9887101
2 -11.107373 0.01 0.0299749 0.8625471
3 -8.901743 0.01 0.0003941 0.9841618
Code 
testdf(variable = financial_data$dy2, #ADF test variable 2
       max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -7.989698 0.01 11.4040782 0.0007328
1 -5.274622 0.01 0.0050530 0.9433307
2 -5.093104 0.01 0.0007863 0.9776301
3 -4.960609 0.01 0.0000187 0.9965522
Code 
       testdf(variable = financial_data$dy3, #ADF test variable 3
              max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -4.529095 0.01 11.9091594 0.0005586
1 -5.539558 0.01 0.0075518 0.9307501
2 -5.387325 0.01 0.0044727 0.9466787
3 -4.529314 0.01 0.0076663 0.9302283
Code 
       testdf(variable = financial_data$dy4, #ADF test variable 4
              max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -14.976767 0.01 0.0087075 0.9256541
1 -10.799791 0.01 0.0003822 0.9844031
2 -8.919976 0.01 0.0067850 0.9343518
3 -5.921353 0.01 0.5225323 0.4697635
Code 
       testdf(variable = financial_data$dy5, #ADF test variable 5
              max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -4.783862 0.01 26.8354772 0.0000002
1 -6.624693 0.01 0.9335859 0.3339329
2 -5.245343 0.01 0.1037083 0.7474241
3 -5.506683 0.01 0.0139857 0.9058609
Code 
       testdf(variable = financial_data$dy6, #ADF test variable 6
              max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -18.012092 0.01 0.0007932 0.9775311
1 -12.542031 0.01 0.0011870 0.9725155
2 -10.204864 0.01 0.0003551 0.9849657
3 -8.516099 0.01 0.0002669 0.9869653
Code 
       testdf(variable = financial_data$dy7, #ADF test variable 7
              max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -8.174416 0.01 12.1264287 0.0004971
1 -5.277486 0.01 0.0030362 0.9560572
2 -5.081266 0.01 0.0004968 0.9822167
3 -5.006448 0.01 0.0003413 0.9852596
Code 
       testdf(variable = financial_data$dy8, #ADF test variable 8
              max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -13.993074 0.01 0.0016383 0.9677134
1 -10.691950 0.01 0.0001454 0.9903780
2 -8.519015 0.01 0.0900690 0.7640895
3 -5.601904 0.01 1.0258157 0.3111435
Code 
       testdf(variable = financial_data$dy9, #ADF test variable 9
              max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -4.763719 0.01 28.6705806 0.0000001
1 -6.682840 0.01 1.1258497 0.2886623
2 -5.224669 0.01 0.1005344 0.7511893
3 -5.511692 0.01 0.0162463 0.8985756
Code 
       testdf(variable = financial_data$dy10, #ADF test variable 10
 
                         max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -6.346090 0.01 6.5970811 0.0102146
1 -4.839809 0.01 0.2597729 0.6102762
2 -5.495769 0.01 0.0529907 0.8179386
3 -4.757224 0.01 0.0030569 0.9559079

P value of Breusch-Godfrey test in 0augmentations is more than significant level(0.05), so there is no autocorrelation and we can check stationary. Morever, p value of ADF test gives us a need for rejecting null hypothesis “unit root”. Our data is ready for regression.

Estimate cointegrate vector

Code 
model.coint <- lm(y3 ~ y10, data = financial_data) #There are 4 cointegrated pairs but I choosed y3 and y10.
summary(model.coint)

Call:
lm(formula = y3 ~ y10, data = financial_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2580 -0.2659 -0.0140  0.2628  1.1026 

Coefficients:
               Estimate  Std. Error t value Pr(>|t|)    
(Intercept) -203.059205    0.241065  -842.3   <2e-16 ***
y10            1.997761    0.001078  1852.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4121 on 298 degrees of freedom
Multiple R-squared:  0.9999,    Adjusted R-squared:  0.9999 
F-statistic: 3.432e+06 on 1 and 298 DF,  p-value: < 2.2e-16

y10 has significant effect on OLS model of Engle Grager method, so one unit increase in y10 increases y3 1.998e+00 unit. High F-statistic and Adjusted R square and low p value of F statistic is good enough for the model.

Checking autocorrelation among residuals

Code 
testdf(variable = residuals(model.coint), max.augmentations = 3)

augmentations adf p_adf bgodfrey p_bg
0 -17.288806 0.01 0.0004325 0.9834088
1 -10.945801 0.01 0.0154520 0.9010729
2 -10.178250 0.01 0.0035422 0.9525409
3 -8.755398 0.01 0.0212039 0.8842248

High negative absolute value of adf in 0 augmentations and lower p-value than 0,05 in ADF test shows that residuals are stationary, so cointegration conditions are satisfied(integration order of residuals is less than integration order of equation itself). Now we will need to esitimate models.

Plot the cointegrated pairs

Code 
plot(financial_data[, c(3,10)],
     col = c("black", "blue"),
     major.ticks = "years", 
     grid.ticks.on = "years",
     grid.ticks.lty = 3,
     main = "",
     legend.loc = "topleft")

Code 
financialdata_y3_y10 <- financial_data[ , c("y3", "y10")] #choosing only cointegrated pairs
tail(financialdata_y3_y10)
               y3    y10
2022-03-22 289.20 246.30
2022-03-23 288.83 245.85
2022-03-24 287.98 245.64
2022-03-25 287.54 245.77
2022-03-26 286.01 244.87
2022-03-27 285.43 244.08

Johansen Cointegration test for find the number of vectors

Code 
johansen_result_trace <- ca.jo(financialdata_y3_y10, type = "trace", K = 6, ecdet = "const")                 
summary(johansen_result_trace) #trace test

###################### 
# Johansen-Procedure # 
###################### 

Test type: trace statistic , without linear trend and constant in cointegration 

Eigenvalues (lambda):
[1] 1.504496e-01 7.861644e-03 7.436528e-18

Values of teststatistic and critical values of test:

          test 10pct  5pct  1pct
r <= 1 |  2.32  7.52  9.24 12.97
r = 0  | 50.26 17.85 19.96 24.60

Eigenvectors, normalised to first column:
(These are the cointegration relations)

             y3.l6     y10.l6   constant
y3.l6      1.00000   1.000000   1.000000
y10.l6    -1.99825  -1.732193  -1.910262
constant 203.17352 140.486475 186.883186

Weights W:
(This is the loading matrix)

          y3.l6       y10.l6               constant
y3.d  0.4125932 -0.013072792 -0.0000000000002899042
y10.d 0.7075202 -0.006356424 -0.0000000000007212011

In first test, test statistic is bigger than 5pct level and we reject Null Hypothesis which shows r=0. Next we test Null Hypothesis r<=1 and cannot reject null. Our conclusion is that there is 1 cointegration vector for pairs. However, we can use alternative of the test.

Alternative variant of the test

Code 
johansen_result_eigen <- ca.jo(financialdata_y3_y10, type = "eigen", K = 6, ecdet = "const") 
summary(johansen_result_eigen) #eigen test

###################### 
# Johansen-Procedure # 
###################### 

Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration 

Eigenvalues (lambda):
[1] 1.504496e-01 7.861644e-03 7.436528e-18

Values of teststatistic and critical values of test:

          test 10pct  5pct  1pct
r <= 1 |  2.32  7.52  9.24 12.97
r = 0  | 47.94 13.75 15.67 20.20

Eigenvectors, normalised to first column:
(These are the cointegration relations)

             y3.l6     y10.l6   constant
y3.l6      1.00000   1.000000   1.000000
y10.l6    -1.99825  -1.732193  -1.910262
constant 203.17352 140.486475 186.883186

Weights W:
(This is the loading matrix)

          y3.l6       y10.l6               constant
y3.d  0.4125932 -0.013072792 -0.0000000000002899042
y10.d 0.7075202 -0.006356424 -0.0000000000007212011
Code 
str(johansen_result_eigen)
Formal class 'ca.jo' [package "urca"] with 26 slots
  ..@ x        : num [1:300, 1:2] 174 174 176 178 181 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:300] "2021-06-01" "2021-06-02" "2021-06-03" "2021-06-04" ...
  .. .. ..$ : chr [1:2] "y3" "y10"
  ..@ Z0       : num [1:294, 1:2] 2.56 1.86 1.51 2.77 2.98 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "y3.d" "y10.d"
  ..@ Z1       : num [1:294, 1:10] 3.34 2.56 1.86 1.51 2.77 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:10] "y3.dl1" "y10.dl1" "y3.dl2" "y10.dl2" ...
  ..@ ZK       : num [1:294, 1:3] 174 174 176 178 181 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:294] "2021-06-01" "2021-06-02" "2021-06-03" "2021-06-04" ...
  .. .. ..$ : chr [1:3] "y3.l6" "y10.l6" "constant"
  ..@ type     : chr "maximal eigenvalue statistic (lambda max)"
  ..@ model    : chr "without linear trend and constant in cointegration"
  ..@ ecdet    : chr "const"
  ..@ lag      : int 6
  ..@ P        : int 2
  ..@ season   : NULL
  ..@ dumvar   : NULL
  ..@ cval     : num [1:2, 1:3] 7.52 13.75 9.24 15.67 12.97 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "r <= 1 |" "r = 0  |"
  .. .. ..$ : chr [1:3] "10pct" "5pct" "1pct"
  ..@ teststat : num [1:2] 2.32 47.94
  ..@ lambda   : num [1:3] 1.50e-01 7.86e-03 7.44e-18
  ..@ Vorg     : num [1:3, 1:3] 5.969 -11.927 1212.676 -0.15 0.259 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:3] "y3.l6" "y10.l6" "constant"
  .. .. ..$ : chr [1:3] "y3.l6" "y10.l6" "constant"
  ..@ V        : num [1:3, 1:3] 1 -2 203.17 1 -1.73 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:3] "y3.l6" "y10.l6" "constant"
  .. .. ..$ : chr [1:3] "y3.l6" "y10.l6" "constant"
  ..@ W        : num [1:2, 1:3] 0.41259316498636 0.70752019008385 -0.01307279161621 -0.00635642359973 -0.00000000000029 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "y3.d" "y10.d"
  .. .. ..$ : chr [1:3] "y3.l6" "y10.l6" "constant"
  ..@ PI       : num [1:2, 1:3] 0.4 0.701 -0.802 -1.403 81.991 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "y3.d" "y10.d"
  .. .. ..$ : chr [1:3] "y3.l6" "y10.l6" "constant"
  ..@ DELTA    : num [1:2, 1:2] 0.988 0.514 0.514 0.307
  ..@ GAMMA    : num [1:2, 1:10] 1.2471 1.1234 -0.3473 -1.179 0.0991 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "y3.d" "y10.d"
  .. .. ..$ : chr [1:10] "y3.dl1" "y10.dl1" "y3.dl2" "y10.dl2" ...
  ..@ R0       : num [1:294, 1:2] -0.2843 -0.0516 -0.2363 1.5016 0.2667 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "R0.y3.d" "R0.y10.d"
  ..@ RK       : num [1:294, 1:3] 142 116 125 119 114 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:294] "2021-06-01" "2021-06-02" "2021-06-03" "2021-06-04" ...
  .. .. ..$ : chr [1:3] "RK.y3.l6" "RK.y10.l6" "RK.constant"
  ..@ bp       : logi NA
  ..@ spec     : chr "longrun"
  ..@ call     : language ca.jo(x = financialdata_y3_y10, type = "eigen", ecdet = "const", K = 6)
  ..@ test.name: chr "Johansen-Procedure"

The alternative variant also gives the same result. There is one cointegrated pair.

Estimatin VECM

Code 
TMS.vec4 <- cajorls(johansen_result_eigen, r = 1) #we can either eigen or trace results
summary(TMS.vec4)
     Length Class  Mode   
rlm  12     mlm    list   
beta  3     -none- numeric
Code 
str(TMS.vec4)
List of 2
 $ rlm :List of 12
  ..$ coefficients : num [1:11, 1:2] 0.413 1.259 -0.355 0.102 -0.577 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:11] "ect1" "y3.dl1" "y10.dl1" "y3.dl2" ...
  .. .. ..$ : chr [1:2] "y3.d" "y10.d"
  ..$ residuals    : num [1:294, 1:2] -0.3069 -0.0262 -0.1902 1.6196 0.3964 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:294] "2021-06-01" "2021-06-02" "2021-06-03" "2021-06-04" ...
  .. .. ..$ : chr [1:2] "y3.d" "y10.d"
  ..$ effects      : num [1:294, 1:2] -2.24 -32.16 -0.18 4.11 0.6 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:294] "ect1" "y3.dl1" "y10.dl1" "y3.dl2" ...
  .. .. ..$ : chr [1:2] "y3.d" "y10.d"
  ..$ rank         : int 11
  ..$ fitted.values: num [1:294, 1:2] 2.87 1.89 1.7 1.15 2.58 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:294] "2021-06-01" "2021-06-02" "2021-06-03" "2021-06-04" ...
  .. .. ..$ : chr [1:2] "y3.d" "y10.d"
  ..$ assign       : int [1:11] 1 2 3 4 5 6 7 8 9 10 ...
  ..$ qr           :List of 5
  .. ..$ qr   : num [1:294, 1:11] -7.0041 0.0476 0.0692 0.0681 -0.1213 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:294] "2021-06-01" "2021-06-02" "2021-06-03" "2021-06-04" ...
  .. .. .. ..$ : chr [1:11] "ect1" "y3.dl1" "y10.dl1" "y3.dl2" ...
  .. .. ..- attr(*, "assign")= int [1:11] 1 2 3 4 5 6 7 8 9 10 ...
  .. ..$ qraux: num [1:11] 1.05 1.06 1.03 1.04 1 ...
  .. ..$ pivot: int [1:11] 1 2 3 4 5 6 7 8 9 10 ...
  .. ..$ tol  : num 0.0000001
  .. ..$ rank : int 11
  .. ..- attr(*, "class")= chr "qr"
  ..$ df.residual  : int 283
  ..$ xlevels      : Named list()
  ..$ call         : language lm(formula = substitute(form1), data = data.mat)
  ..$ terms        :Classes 'terms', 'formula'  language z@Z0 ~ ect1 + y3.dl1 + y10.dl1 + y3.dl2 + y10.dl2 + y3.dl3 + y10.dl3 +      y3.dl4 + y10.dl4 + y3.dl5 + y10.dl5 - 1
  .. .. ..- attr(*, "variables")= language list(z@Z0, ect1, y3.dl1, y10.dl1, y3.dl2, y10.dl2, y3.dl3, y10.dl3, y3.dl4,      y10.dl4, y3.dl5, y10.dl5)
  .. .. ..- attr(*, "factors")= int [1:12, 1:11] 0 1 0 0 0 0 0 0 0 0 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:12] "z@Z0" "ect1" "y3.dl1" "y10.dl1" ...
  .. .. .. .. ..$ : chr [1:11] "ect1" "y3.dl1" "y10.dl1" "y3.dl2" ...
  .. .. ..- attr(*, "term.labels")= chr [1:11] "ect1" "y3.dl1" "y10.dl1" "y3.dl2" ...
  .. .. ..- attr(*, "order")= int [1:11] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. ..- attr(*, "intercept")= int 0
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: 0x00000242b1b729a8> 
  .. .. ..- attr(*, "predvars")= language list(z@Z0, ect1, y3.dl1, y10.dl1, y3.dl2, y10.dl2, y3.dl3, y10.dl3, y3.dl4,      y10.dl4, y3.dl5, y10.dl5)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:12] "nmatrix.2" "numeric" "numeric" "numeric" ...
  .. .. .. ..- attr(*, "names")= chr [1:12] "z@Z0" "ect1" "y3.dl1" "y10.dl1" ...
  ..$ model        :'data.frame':   294 obs. of  12 variables:
  .. ..$ z@Z0   : num [1:294, 1:2] 2.56 1.86 1.51 2.77 2.98 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : NULL
  .. .. .. ..$ : chr [1:2] "y3.d" "y10.d"
  .. ..$ ect1   : num [1:294] 0.373 0.333 0.485 0.477 -0.849 ...
  .. ..$ y3.dl1 : num [1:294] 3.34 2.56 1.86 1.51 2.77 ...
  .. ..$ y10.dl1: num [1:294] 1.41 1.3 0.77 0.96 0.99 ...
  .. ..$ y3.dl2 : num [1:294] 3.21 3.34 2.56 1.86 1.51 ...
  .. ..$ y10.dl2: num [1:294] 2.27 1.41 1.3 0.77 0.96 ...
  .. ..$ y3.dl3 : num [1:294] 2.19 3.21 3.34 2.56 1.86 ...
  .. ..$ y10.dl3: num [1:294] 1.1 2.27 1.41 1.3 0.77 ...
  .. ..$ y3.dl4 : num [1:294] 2.13 2.19 3.21 3.34 2.56 ...
  .. ..$ y10.dl4: num [1:294] 0.99 1.1 2.27 1.41 1.3 ...
  .. ..$ y3.dl5 : num [1:294] 0 2.13 2.19 3.21 3.34 ...
  .. ..$ y10.dl5: num [1:294] 0.02 0.99 1.1 2.27 1.41 ...
  .. ..- attr(*, "terms")=Classes 'terms', 'formula'  language z@Z0 ~ ect1 + y3.dl1 + y10.dl1 + y3.dl2 + y10.dl2 + y3.dl3 + y10.dl3 +      y3.dl4 + y10.dl4 + y3.dl5 + y10.dl5 - 1
  .. .. .. ..- attr(*, "variables")= language list(z@Z0, ect1, y3.dl1, y10.dl1, y3.dl2, y10.dl2, y3.dl3, y10.dl3, y3.dl4,      y10.dl4, y3.dl5, y10.dl5)
  .. .. .. ..- attr(*, "factors")= int [1:12, 1:11] 0 1 0 0 0 0 0 0 0 0 ...
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:12] "z@Z0" "ect1" "y3.dl1" "y10.dl1" ...
  .. .. .. .. .. ..$ : chr [1:11] "ect1" "y3.dl1" "y10.dl1" "y3.dl2" ...
  .. .. .. ..- attr(*, "term.labels")= chr [1:11] "ect1" "y3.dl1" "y10.dl1" "y3.dl2" ...
  .. .. .. ..- attr(*, "order")= int [1:11] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. .. ..- attr(*, "intercept")= int 0
  .. .. .. ..- attr(*, "response")= int 1
  .. .. .. ..- attr(*, ".Environment")=<environment: 0x00000242b1b729a8> 
  .. .. .. ..- attr(*, "predvars")= language list(z@Z0, ect1, y3.dl1, y10.dl1, y3.dl2, y10.dl2, y3.dl3, y10.dl3, y3.dl4,      y10.dl4, y3.dl5, y10.dl5)
  .. .. .. ..- attr(*, "dataClasses")= Named chr [1:12] "nmatrix.2" "numeric" "numeric" "numeric" ...
  .. .. .. .. ..- attr(*, "names")= chr [1:12] "z@Z0" "ect1" "y3.dl1" "y10.dl1" ...
  ..- attr(*, "class")= chr [1:2] "mlm" "lm"
 $ beta: num [1:3, 1] 1 -2 203
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "y3.l6" "y10.l6" "constant"
  .. ..$ : chr "ect1"
Code 
summary(TMS.vec4$rlm)
Response y3.d :

Call:
lm(formula = y3.d ~ ect1 + y3.dl1 + y10.dl1 + y3.dl2 + y10.dl2 + 
    y3.dl3 + y10.dl3 + y3.dl4 + y10.dl4 + y3.dl5 + y10.dl5 - 
    1, data = data.mat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.69689 -0.61031  0.07821  0.64112  2.93975 

Coefficients:
        Estimate Std. Error t value          Pr(>|t|)    
ect1     0.41259    0.35408   1.165             0.245    
y3.dl1   1.25943    0.16609   7.583 0.000000000000488 ***
y10.dl1 -0.35514    0.29790  -1.192             0.234    
y3.dl2   0.10212    0.22958   0.445             0.657    
y10.dl2 -0.57742    0.41891  -1.378             0.169    
y3.dl3   0.04274    0.26299   0.163             0.871    
y10.dl3 -0.45330    0.49456  -0.917             0.360    
y3.dl4   0.24502    0.30164   0.812             0.417    
y10.dl4 -0.11653    0.57378  -0.203             0.839    
y3.dl5   0.23621    0.32386   0.729             0.466    
y10.dl5 -0.56602    0.64404  -0.879             0.380    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.017 on 283 degrees of freedom
Multiple R-squared:  0.7848,    Adjusted R-squared:  0.7764 
F-statistic: 93.82 on 11 and 283 DF,  p-value: < 2.2e-16


Response y10.d :

Call:
lm(formula = y10.d ~ ect1 + y3.dl1 + y10.dl1 + y3.dl2 + y10.dl2 + 
    y3.dl3 + y10.dl3 + y3.dl4 + y10.dl4 + y3.dl5 + y10.dl5 - 
    1, data = data.mat)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.72167 -0.34499 -0.00932  0.38694  1.60174 

Coefficients:
        Estimate Std. Error t value         Pr(>|t|)    
ect1     0.70752    0.19705   3.591         0.000389 ***
y3.dl1   1.12945    0.09243  12.220          < 2e-16 ***
y10.dl1 -1.18283    0.16578  -7.135 0.00000000000814 ***
y3.dl2   0.49607    0.12776   3.883         0.000129 ***
y10.dl2 -1.18484    0.23312  -5.083 0.00000067784059 ***
y3.dl3   0.52732    0.14635   3.603         0.000371 ***
y10.dl3 -1.21963    0.27522  -4.431 0.00001339675361 ***
y3.dl4   0.57103    0.16786   3.402         0.000766 ***
y10.dl4 -1.01878    0.31930  -3.191         0.001579 ** 
y3.dl5   0.62640    0.18023   3.476         0.000590 ***
y10.dl5 -1.24639    0.35840  -3.478         0.000585 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5661 on 283 degrees of freedom
Multiple R-squared:  0.7536,    Adjusted R-squared:  0.744 
F-statistic: 78.67 on 11 and 283 DF,  p-value: < 2.2e-16
Code 
TMS.vec4$beta # Cointegarting vector
              ect1
y3.l6      1.00000
y10.l6    -1.99825
constant 203.17352

There 5 cointegration matrices. P values of Fstatistic is less 0.05, so we reject null hypothesis that coefficents are 0. Morever, ect1 in y10.d is 0.70752 and shows p value<0.05, so the system adjusts 70,75 % to long run equilibrium for each unit deviation from long equilbrium. On the other hand, because ect1 is not statistically significant in y3.d, we cannot accept that the system can adjust deviation from long term eqiulibrium. $beta code gave use cointegration vector.

Repramatrize from VECM to VAR

Code 
TMS.vec4.asVAR <- vec2var(johansen_result_eigen, r = 1)
TMS.vec4.asVAR

Coefficient matrix of lagged endogenous variables:

A1:
       y3.l1     y10.l1
y3  2.259433 -0.3551365
y10 1.129454 -0.1828322


A2:
         y3.l2       y10.l2
y3  -1.1573099 -0.222285118
y10 -0.6333811 -0.002007442


A3:
          y3.l3      y10.l3
y3  -0.05938044  0.12412110
y10  0.03124713 -0.03478676


A4:
         y3.l4    y10.l4
y3  0.20227998 0.3367707
y10 0.04370972 0.2008484


A5:
           y3.l5     y10.l5
y3  -0.008811775 -0.4494859
y10  0.055367291 -0.2276139


A6:
         y3.l6     y10.l6
y3  0.17638221 -0.2584488
y10 0.08112279 -0.1674107


Coefficient matrix of deterministic regressor(s).

     constant
y3   83.82801
y10 143.74937

Because generally VECM is differenced equation, VAR has one more matrices. So VAR(p+1)=VECM. Here we have 6 matrices.

VAR diagnostics in IRF

Code 
irf_result <- irf(TMS.vec4.asVAR)
plot(irf_result)

We can see that shocks coming from y3 has more impact on both and after 10 periods these shocks disappear. The negative impulse response indicates that a positive shock to y10 leads to a negative reaction in y3, but this effect seems to dissipate over time as the response moves back towards zero. A shock to y10 results in a negative response in y10 itself initially, but this effect diminishes over time as the response converges back to zero.

VAR diagnostics on FEVD

Code 
plot(fevd(TMS.vec4.asVAR, n.ahead = 20))

Shocks to y10 account for nearly all the forecast error variance, indicating a strong dependence of y3 on y10. After the first period, y10 almost entirely explains its own forecast error variance, highlighting its self-determination in forecasting errors.

Perform the Portmanteau test for residuals on VAR model

Code 
head(residuals(TMS.vec4.asVAR))
     resids of y3 resids of y10
[1,]  -0.30693742    0.03265190
[2,]  -0.02622717    0.04001272
[3,]  -0.19019950    0.13458942
[4,]   1.61957358    0.65575761
[5,]   0.39639475    0.42779874
[6,]   1.06454554    0.55720276
Code 
portmanteau_test <- serial.test(TMS.vec4.asVAR)
print(portmanteau_test)

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object TMS.vec4.asVAR
Chi-squared = 40.023, df = 42, p-value = 0.5581
$serial

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object TMS.vec4.asVAR
Chi-squared = 40.023, df = 42, p-value = 0.5581

We cannot reject Null Hypothesis which indicates no autororeelation in residuals, so there is no autocorrelation among residuals.

Perform the Jarque-Bera test for normality

Code 
jb_test <- normality.test(TMS.vec4.asVAR)
print(jb_test)
$JB

    JB-Test (multivariate)

data:  Residuals of VAR object TMS.vec4.asVAR
Chi-squared = 0.62506, df = 4, p-value = 0.9602


$Skewness

    Skewness only (multivariate)

data:  Residuals of VAR object TMS.vec4.asVAR
Chi-squared = 0.24287, df = 2, p-value = 0.8857


$Kurtosis

    Kurtosis only (multivariate)

data:  Residuals of VAR object TMS.vec4.asVAR
Chi-squared = 0.38219, df = 2, p-value = 0.8261
$jb.mul
$jb.mul$JB

    JB-Test (multivariate)

data:  Residuals of VAR object TMS.vec4.asVAR
Chi-squared = 0.62506, df = 4, p-value = 0.9602


$jb.mul$Skewness

    Skewness only (multivariate)

data:  Residuals of VAR object TMS.vec4.asVAR
Chi-squared = 0.24287, df = 2, p-value = 0.8857


$jb.mul$Kurtosis

    Kurtosis only (multivariate)

data:  Residuals of VAR object TMS.vec4.asVAR
Chi-squared = 0.38219, df = 2, p-value = 0.8261

P values are bigger than 0.05, so we cannot ignore Null Hypothesis which shows normal distribution. Our residuals are normally distributed.

Forecasting out of sample

Code 
TMS.testdata <- financial_data["/3/7/2022", ]
TMS.testdata <- subset(TMS.testdata, select = c(y3, y10))
tail(TMS.testdata)
               y3    y10
2022-03-22 289.20 246.30
2022-03-23 288.83 245.85
2022-03-24 287.98 245.64
2022-03-25 287.54 245.77
2022-03-26 286.01 244.87
2022-03-27 285.43 244.08
Code 
johan.test.eigen2.testdata <- ca.jo(TMS.testdata,
                                 ecdet = "const",
                                 type = "trace", 
                                 K = 4)
summary(johan.test.eigen2.testdata)

###################### 
# Johansen-Procedure # 
###################### 

Test type: trace statistic , without linear trend and constant in cointegration 

Eigenvalues (lambda):
[1] 2.126143e-01 9.370866e-03 1.734723e-18

Values of teststatistic and critical values of test:

          test 10pct  5pct  1pct
r <= 1 |  2.79  7.52  9.24 12.97
r = 0  | 73.54 17.85 19.96 24.60

Eigenvectors, normalised to first column:
(These are the cointegration relations)

              y3.l4     y10.l4  constant
y3.l4      1.000000   1.000000  1.000000
y10.l4    -1.998201  -1.724705 -1.304082
constant 203.166419 138.259152 71.888373

Weights W:
(This is the loading matrix)

          y3.l4       y10.l4      constant
y3.d  0.1260953 -0.013703907 -8.908024e-16
y10.d 0.5689688 -0.006982952 -5.897452e-16

For out of sample, we tested Johansen test again and repramatrized it. Here there is only one vector.

Code 
TMS.vec4.fore <- 
  predict(
    vec2var(
      johan.test.eigen2.testdata, 
      r = 1),     # no of cointegrating vectors 
    n.ahead = 20, # forecast horizon
    ci = 0.95)    # confidence level for intervals

TMS.vec4.fore$fcst$y3
          fcst    lower    upper        CI
 [1,] 285.2120 283.2259 287.1980  1.986034
 [2,] 285.2980 280.7385 289.8576  4.559572
 [3,] 285.4969 277.9828 293.0110  7.514084
 [4,] 285.5471 274.9346 296.1595 10.612434
 [5,] 285.5513 271.8319 299.2706 13.719334
 [6,] 285.5234 268.7630 302.2838 16.760351
 [7,] 285.5071 265.8111 305.2031 19.695973
 [8,] 285.4938 262.9869 308.0007 22.506879
 [9,] 285.4864 260.2999 310.6730 25.186535
[10,] 285.4796 257.7435 313.2158 27.736139
[11,] 285.4746 255.3132 315.6359 30.161362
[12,] 285.4703 253.0002 317.9403 32.470075
[13,] 285.4671 250.7961 320.1381 34.671019
[14,] 285.4647 248.6917 322.2377 36.773010
[15,] 285.4629 246.6784 324.2475 38.784517
[16,] 285.4616 244.7481 326.1750 40.713444
[17,] 285.4606 242.8936 328.0276 42.567036
[18,] 285.4598 241.1080 329.8117 44.351862
[19,] 285.4593 239.3854 331.5331 46.073827
[20,] 285.4589 237.7206 333.1971 47.738225
Code 
TMS.vec4.fore$fcst$y10
          fcst    lower    upper        CI
 [1,] 244.4091 243.3038 245.5143  1.105265
 [2,] 244.4101 242.0766 246.7435  2.333497
 [3,] 244.5979 240.8117 248.3840  3.786132
 [4,] 244.5708 239.2513 249.8904  5.319554
 [5,] 244.5863 237.7225 251.4502  6.863865
 [6,] 244.5589 236.1792 252.9386  8.379684
 [7,] 244.5581 234.7125 254.4038  9.845634
 [8,] 244.5489 233.2982 255.7995 11.250646
 [9,] 244.5473 231.9564 257.1382 12.590880
[10,] 244.5428 230.6763 258.4093 13.866521
[11,] 244.5407 229.4605 259.6210 15.080244
[12,] 244.5382 228.3024 260.7741 16.235844
[13,] 244.5368 227.1992 261.8744 17.337624
[14,] 244.5355 226.1456 262.9254 18.389937
[15,] 244.5347 225.1377 263.9317 19.396991
[16,] 244.5340 224.1712 264.8967 20.362721
[17,] 244.5335 223.2427 265.8242 21.290742
[18,] 244.5331 222.3488 266.7174 22.184333
[19,] 244.5328 221.4864 267.5793 23.046447
[20,] 244.5326 220.6529 268.4123 23.879731
Code 
tail(index(financial_data), 20)
 [1] "2022-03-08" "2022-03-09" "2022-03-10" "2022-03-11" "2022-03-12"
 [6] "2022-03-13" "2022-03-14" "2022-03-15" "2022-03-16" "2022-03-17"
[11] "2022-03-18" "2022-03-19" "2022-03-20" "2022-03-21" "2022-03-22"
[16] "2022-03-23" "2022-03-24" "2022-03-25" "2022-03-26" "2022-03-27"

Here we added forecast series and lower and upper bounds.

Adding y3 and y10 forecasts

Code 
y3_forecast <- xts(TMS.vec4.fore$fcst$y3[,-4], 
                   # we exclude the last column with CI
                   tail(index(financial_data), 20))

names(y3_forecast) <- c("y3_fore", "y3_lower", "y3_upper")

y10_forecast <- xts(TMS.vec4.fore$fcst$y10[,-4],
                   # we exclude the last column with CI
                   tail(index(financial_data), 20))
names(y10_forecast) <- c("y10_fore", "y10_lower", "y10_upper")

Combinening real data and forecast data

Code 
financial_data_fore <- merge(financial_data, 
                 y3_forecast,
                 y10_forecast)

Plotting forecats

Code 
plot(financial_data_fore[index(financial_data_fore) > as.Date("2022-03-07"), c("y3", "y3_fore",
                                             "y3_lower", "y3_upper")], 
     
     main = "20 days forecast of instrument y3",
     col = c("black", "blue", "red", "red"))

Code 
plot(financial_data_fore[index(financial_data_fore) > as.Date("2022-03-07"), c("y10", "y10_fore",
                                                                     "y10_lower", "y10_upper")], 
     
     main = "20 days forecast of instrument y10",
     col = c("black", "blue", "red", "red"))

The forecasted prices for y3 and y10 over 20 time periods show that projected observations fall within the upper and lower boundaries. In the plot, the blue line represents the forecasted values, while the other line (assumed to be black) indicates the original observations.

Calculating forecast errors for y3

Code 
financial_data_fore$mae.y3   <-  abs(financial_data_fore$y3 - financial_data_fore$y3_fore)
tail(financial_data_fore$mae.y3, 20)
                mae.y3
2022-03-08 17.67801884
2022-03-09 14.81196569
2022-03-10 12.72307274
2022-03-11 11.56294243
2022-03-12 11.10873033
2022-03-13 11.30660064
2022-03-14 10.53288456
2022-03-15  8.32617952
2022-03-16  6.00355962
2022-03-17  4.57037360
2022-03-18  3.32543082
2022-03-19  2.90974191
2022-03-20  2.97287567
2022-03-21  3.39530022
2022-03-22  3.73706357
2022-03-23  3.36840866
2022-03-24  2.51940563
2022-03-25  2.08016109
2022-03-26  0.55072322
2022-03-27  0.02885409
Code 
financial_data_fore$mse.y3   <-  (financial_data_fore$y3 - financial_data_fore$y3_fore) ^ 2
tail(financial_data_fore$mse.y3, 20)
                   mse.y3
2022-03-08 312.5123502194
2022-03-09 219.3943275444
2022-03-10 161.8765799741
2022-03-11 133.7016375905
2022-03-12 123.4038895259
2022-03-13 127.8392180972
2022-03-14 110.9416571973
2022-03-15  69.3252654070
2022-03-16  36.0427280663
2022-03-17  20.8883148576
2022-03-18  11.0584901296
2022-03-19   8.4665979558
2022-03-20   8.8379897690
2022-03-21  11.5280635822
2022-03-22  13.9656440971
2022-03-23  11.3461768767
2022-03-24   6.3474047119
2022-03-25   4.3270701417
2022-03-26   0.3032960660
2022-03-27   0.0008325584
Code 
financial_data_fore$mape.y3   <-  abs((financial_data_fore$y3 - financial_data_fore$y3_fore) / financial_data_fore$y3_fore)
tail(financial_data_fore$amape.y3, 20)
NULL
Code 
financial_data_fore$amape.y3   <-  abs((financial_data_fore$y3 - financial_data_fore$y3_fore) / (financial_data_fore$y3+financial_data_fore$y3_fore))
tail(financial_data_fore$amape.y3, 20)
                amape.y3
2022-03-08 0.03005944447
2022-03-09 0.02530195149
2022-03-10 0.02179664859
2022-03-11 0.01984519415
2022-03-12 0.01908023927
2022-03-13 0.01941535957
2022-03-14 0.01811183356
2022-03-15 0.01437248556
2022-03-16 0.01040520755
2022-03-17 0.00794116131
2022-03-18 0.00579066356
2022-03-19 0.00507055955
2022-03-20 0.00518006407
2022-03-21 0.00591181299
2022-03-22 0.00650305306
2022-03-23 0.00586532818
2022-03-24 0.00439349019
2022-03-25 0.00363029960
2022-03-26 0.00096369699
2022-03-27 0.00005054239

Calculate means for the last 20 values for each error metric and store them in a table

Code 
error_metrics_means <- data.frame(
  Metric = c("amape.y3", "mape.y3", "mse.y3", "mae.y3"),
  Mean = c(
    mean(tail(financial_data_fore$amape.y3, 20), na.rm = TRUE),
    mean(tail(financial_data_fore$mape.y3, 20), na.rm = TRUE),
    mean(tail(financial_data_fore$mse.y3, 20), na.rm = TRUE),
    mean(tail(financial_data_fore$mae.y3, 20), na.rm = TRUE)
  )
)

Calculating forecast errors for y10

Code 
financial_data_fore$mae.y10   <-  abs(financial_data_fore$y10 - financial_data_fore$y10_fore)
tail(financial_data_fore$mae.y10, 20)
             mae.y10
2022-03-08 8.9209393
2022-03-09 7.1899471
2022-03-10 6.2421233
2022-03-11 5.5691888
2022-03-12 5.3336515
2022-03-13 5.7511118
2022-03-14 5.1918711
2022-03-15 4.4911458
2022-03-16 2.5427150
2022-03-17 2.1972097
2022-03-18 1.6892765
2022-03-19 1.8817754
2022-03-20 1.8431908
2022-03-21 2.1444948
2022-03-22 1.7653333
2022-03-23 1.3160375
2022-03-24 1.1065269
2022-03-25 1.2369159
2022-03-26 0.3371955
2022-03-27 0.4525891
Code 
financial_data_fore$mse.y10   <-  (financial_data_fore$y10 - financial_data_fore$y10_fore) ^ 2
tail(financial_data_fore$mse.y10, 20)
              mse.y10
2022-03-08 79.5831578
2022-03-09 51.6953393
2022-03-10 38.9641028
2022-03-11 31.0158640
2022-03-12 28.4478382
2022-03-13 33.0752864
2022-03-14 26.9555257
2022-03-15 20.1703910
2022-03-16  6.4653994
2022-03-17  4.8277305
2022-03-18  2.8536552
2022-03-19  3.5410788
2022-03-20  3.3973524
2022-03-21  4.5988577
2022-03-22  3.1164017
2022-03-23  1.7319546
2022-03-24  1.2244018
2022-03-25  1.5299610
2022-03-26  0.1137008
2022-03-27  0.2048369
Code 
financial_data_fore$mape.y10   <-  abs((financial_data_fore$y10 - financial_data_fore$y10_fore) / financial_data_fore$y10_fore)
tail(financial_data_fore$mape.y10, 20)
              mape.y10
2022-03-08 0.036500035
2022-03-09 0.029417559
2022-03-10 0.025519941
2022-03-11 0.022771273
2022-03-12 0.021806824
2022-03-13 0.023516266
2022-03-14 0.021229599
2022-03-15 0.018365025
2022-03-16 0.010397641
2022-03-17 0.008984970
2022-03-18 0.006907956
2022-03-19 0.007695220
2022-03-20 0.007537478
2022-03-21 0.008769666
2022-03-22 0.007219154
2022-03-23 0.005381819
2022-03-24 0.004525053
2022-03-25 0.005058276
2022-03-26 0.001378938
2022-03-27 0.001850833
Code 
financial_data_fore$amape.y10   <-  abs((financial_data_fore$y10 - financial_data_fore$y10_fore) / (financial_data_fore$y10+financial_data_fore$y10_fore))
tail(financial_data_fore$amape.y10, 20)
              amape.y10
2022-03-08 0.0179229239
2022-03-09 0.0144955673
2022-03-10 0.0125992048
2022-03-11 0.0112574633
2022-03-12 0.0107858099
2022-03-13 0.0116214858
2022-03-14 0.0105033092
2022-03-15 0.0090989612
2022-03-16 0.0051719327
2022-03-17 0.0044723930
2022-03-18 0.0034420890
2022-03-19 0.0038328626
2022-03-20 0.0037545890
2022-03-21 0.0043656903
2022-03-22 0.0035965946
2022-03-23 0.0026836878
2022-03-24 0.0022574190
2022-03-25 0.0025227578
2022-03-26 0.0006889937
2022-03-27 0.0009262739

Calculate means for the last 20 values for each error metric and store them in a table

Code 
error_metrics_means <- data.frame(
  Metric = c("amape.y10", "mape.y10", "mse.y10", "mae.y10"),
  Mean = c(
    mean(tail(financial_data_fore$amape.y10, 20), na.rm = TRUE),
    mean(tail(financial_data_fore$mape.y10, 20), na.rm = TRUE),
    mean(tail(financial_data_fore$mse.y10, 20), na.rm = TRUE),
    mean(tail(financial_data_fore$mae.y10, 20), na.rm = TRUE)
  )
)

Arima modelling

Code 
head(financial_data)
               y1     y2     y3     y4     y5     y6     y7     y8     y9
2021-06-01 176.04 173.57 173.57 188.66 183.66 178.78 168.46 173.57 173.57
2021-06-02 173.37 173.57 173.57 188.46 183.52 181.10 168.41 173.57 173.57
2021-06-03 173.42 173.57 175.70 188.57 183.75 187.25 168.47 173.57 173.57
2021-06-04 178.85 174.52 177.89 188.78 182.79 190.82 170.08 173.57 173.18
2021-06-05 175.39 174.46 181.10 188.11 185.46 187.30 170.15 172.94 174.50
2021-06-06 173.89 174.54 184.44 187.61 190.18 195.26 170.33 171.78 176.99
              y10   dy1   dy2  dy3   dy4   dy5   dy6   dy7   dy8   dy9 dy10
2021-06-01 188.35    NA    NA   NA    NA    NA    NA    NA    NA    NA   NA
2021-06-02 188.37 -2.67  0.00 0.00 -0.20 -0.14  2.32 -0.05  0.00  0.00 0.02
2021-06-03 189.36  0.05  0.00 2.13  0.11  0.23  6.15  0.06  0.00  0.00 0.99
2021-06-04 190.46  5.43  0.95 2.19  0.21 -0.96  3.57  1.61  0.00 -0.39 1.10
2021-06-05 192.73 -3.46 -0.06 3.21 -0.67  2.67 -3.52  0.07 -0.63  1.32 2.27
2021-06-06 194.14 -1.50  0.08 3.34 -0.50  4.72  7.96  0.18 -1.16  2.49 1.41
Code 
str(financial_data)
An xts object on 2021-06-01 / 2022-03-27 containing: 
  Data:    double [300, 20]
  Columns: y1, y2, y3, y4, y5 ... with 15 more columns
  Index:   Date [300] (TZ: "UTC")
Code 
tail(financial_data)
               y1     y2     y3     y4    y5    y6     y7     y8     y9    y10
2022-03-22 268.68 196.26 289.20 174.51 36.93 97.94 204.83 145.53 100.31 246.30
2022-03-23 266.55 196.34 288.83 174.58 36.87 87.23 205.03 145.76 100.22 245.85
2022-03-24 266.38 196.25 287.98 174.18 35.33 78.05 205.08 144.81  99.39 245.64
2022-03-25 272.44 197.49 287.54 174.47 34.93 85.77 206.84 144.90  99.22 245.77
2022-03-26 267.64 198.03 286.01 173.46 35.01 95.43 208.08 143.75  99.22 244.87
2022-03-27 264.95 198.61 285.43 174.40 33.51 95.25 208.76 145.46  98.49 244.08
             dy1   dy2   dy3   dy4   dy5    dy6   dy7   dy8   dy9  dy10
2022-03-22 -0.78 -0.30  0.34 -0.94  0.22   0.84 -0.56 -2.26  0.14 -0.38
2022-03-23 -2.13  0.08 -0.37  0.07 -0.06 -10.71  0.20  0.23 -0.09 -0.45
2022-03-24 -0.17 -0.09 -0.85 -0.40 -1.54  -9.18  0.05 -0.95 -0.83 -0.21
2022-03-25  6.06  1.24 -0.44  0.29 -0.40   7.72  1.76  0.09 -0.17  0.13
2022-03-26 -4.80  0.54 -1.53 -1.01  0.08   9.66  1.24 -1.15  0.00 -0.90
2022-03-27 -2.69  0.58 -0.58  0.94 -1.50  -0.18  0.68  1.71 -0.73 -0.79
Code 
TMS.arimapart <- subset(financial_data, select = c(y3, y10, dy3, dy10))
plot(TMS.arimapart, 
     # plot data in two panels (each column separately)
     multi.panel = 2,
     main = "Original and differenced data",
     col = c("darkgreen", "darkblue"),
     major.ticks = "years", 
     grid.ticks.on = "years",
     grid.ticks.lty = 3,
     yaxis.same = FALSE,
     cex = .6, lwd = 1)

Box Jenkins Procedure for y3

Code 
par(mfrow = c(2, 1)) 
acf(TMS.arimapart$dy3,
    lag.max = 36, # max lag for ACF
    ylim = c(-0.5, 0.5),   # limits for the y axis - we give c(min, max)
    lwd = 5,               # line width
    col = "dark green",
    na.action = na.pass)   # do not stop if there are missing values in the data
pacf(TMS.arimapart$dy3, 
     lag.max = 36, 
     ylim = c(-0.5, 0.5),   # limits for the y axis - we give c(min, max)
     lwd = 5, col = "dark green",
     na.action = na.pass)

For y3, ACF values are significant till 10th lag and PACF values are significant till 2.

Box Jenkins procedure for y10

Code 
par(mfrow = c(2, 1)) 
acf(TMS.arimapart$dy10,
    lag.max = 36, # max lag for ACF
    ylim = c(-0.5, 0.5),   # limits for the y axis - we give c(min, max)
    lwd = 5,               # line width
    col = "dark green",
    na.action = na.pass)   # do not stop if there are missing values in the data
pacf(TMS.arimapart$dy10, 
     lag.max = 36, 
     lwd = 5, col = "dark green",
     na.action = na.pass)

Code 
par(mfrow = c(1, 1))

For y10, ACF values till 10th lag and PACF values till 3rd lag are significant.

Checking models1 for y3

Code 
arima111 <- Arima(TMS.arimapart$y3,  # variable
                  order = c(1, 1, 1)  # (p,d,q) parameters
)
arima111
Series: TMS.arimapart$y3 
ARIMA(1,1,1) 

Coefficients:
         ar1     ma1
      0.8232  0.2170
s.e.  0.0356  0.0558

sigma^2 = 1.059:  log likelihood = -432.55
AIC=871.09   AICc=871.17   BIC=882.19
Code 
coeftest(arima111) #shows whether coefficenta are signifant or not

z test of coefficients:

    Estimate Std. Error z value   Pr(>|z|)    
ar1 0.823240   0.035604 23.1224  < 2.2e-16 ***
ma1 0.217016   0.055754  3.8924 0.00009926 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
summary(arima111) 
Series: TMS.arimapart$y3 
ARIMA(1,1,1) 

Coefficients:
         ar1     ma1
      0.8232  0.2170
s.e.  0.0356  0.0558

sigma^2 = 1.059:  log likelihood = -432.55
AIC=871.09   AICc=871.17   BIC=882.19

Training set error measures:
                     ME     RMSE       MAE        MPE      MAPE     MASE
Training set 0.05390838 1.023838 0.8072486 0.02915556 0.3497714 0.457447
                   ACF1
Training set 0.02136456
Code 
arima111_2 <- Arima(TMS.arimapart$y3,  # variable
                    order = c(1, 1, 1),  # (p,d,q) parameters
                    include.constant = TRUE)  # including a constant
arima111_2
Series: TMS.arimapart$y3 
ARIMA(1,1,1) with drift 

Coefficients:
         ar1     ma1   drift
      0.8179  0.2188  0.3510
s.e.  0.0362  0.0558  0.3903

sigma^2 = 1.06:  log likelihood = -432.15
AIC=872.31   AICc=872.44   BIC=887.11
Code 
coeftest(arima111_2)

z test of coefficients:

      Estimate Std. Error z value   Pr(>|z|)    
ar1   0.817850   0.036232 22.5728  < 2.2e-16 ***
ma1   0.218839   0.055793  3.9224 0.00008768 ***
drift 0.351044   0.390322  0.8994     0.3685    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
#drift is not significant we ignore it.

Checking models2 for y3

Code 
arima112 <- Arima(TMS.arimapart$y3, # variable
                   order = c(1,1,2 ))
arima112
Series: TMS.arimapart$y3 
ARIMA(1,1,2) 

Coefficients:
         ar1     ma1     ma2
      0.7421  0.3516  0.1970
s.e.  0.0568  0.0832  0.0766

sigma^2 = 1.04:  log likelihood = -429.41
AIC=866.82   AICc=866.96   BIC=881.62
Code 
coeftest(arima112)

z test of coefficients:

    Estimate Std. Error z value   Pr(>|z|)    
ar1 0.742066   0.056838 13.0557  < 2.2e-16 ***
ma1 0.351641   0.083159  4.2285 0.00002352 ***
ma2 0.197019   0.076590  2.5724     0.0101 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
arima112_2 <- Arima(TMS.arimapart$y3,  # variable
                    order = c(1, 1, 2),  # (p,d,q) parameters
                    include.constant = TRUE)  # including a constant
arima112_2
Series: TMS.arimapart$y3 
ARIMA(1,1,2) with drift 

Coefficients:
         ar1     ma1     ma2   drift
      0.7324  0.3579  0.2019  0.3541
s.e.  0.0583  0.0837  0.0767  0.3379

sigma^2 = 1.04:  log likelihood = -428.88
AIC=867.76   AICc=867.96   BIC=886.26
Code 
coeftest(arima112_2)

z test of coefficients:

      Estimate Std. Error z value   Pr(>|z|)    
ar1   0.732442   0.058279 12.5678  < 2.2e-16 ***
ma1   0.357874   0.083705  4.2754 0.00001908 ***
ma2   0.201922   0.076706  2.6324   0.008478 ** 
drift 0.354075   0.337886  1.0479   0.294678    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
#drift is not significant we ignore it.

It has lower information criteria and all of coefficents are significant. At least, it is better than ARIMA(1,1,1).

Checking models3 for y3

Code 
arima211 <- Arima(TMS.arimapart$y3,  # variable
                  order = c(2, 1, 1)  # (p,d,q) parameters
)
arima211
Series: TMS.arimapart$y3 
ARIMA(2,1,1) 

Coefficients:
         ar1      ar2     ma1
      1.1163  -0.2679  -0.047
s.e.  0.1744   0.1540   0.176

sigma^2 = 1.053:  log likelihood = -431.3
AIC=870.6   AICc=870.73   BIC=885.4
Code 
coeftest(arima211)

z test of coefficients:

     Estimate Std. Error z value        Pr(>|z|)    
ar1  1.116331   0.174355  6.4026 0.0000000001527 ***
ar2 -0.267948   0.153973 -1.7402         0.08182 .  
ma1 -0.047011   0.176040 -0.2670         0.78943    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
arima211_2 <- Arima(TMS.arimapart$y3,  # variable
                    order = c(2, 1, 1),  # (p,d,q) parameters
                    include.constant = TRUE)  # including a constant
arima211_2
Series: TMS.arimapart$y3 
ARIMA(2,1,1) with drift 

Coefficients:
         ar1      ar2      ma1   drift
      1.1207  -0.2762  -0.0542  0.3523
s.e.  0.1738   0.1529   0.1760  0.3545

sigma^2 = 1.054:  log likelihood = -430.82
AIC=871.64   AICc=871.85   BIC=890.15
Code 
coeftest(arima211_2)

z test of coefficients:

       Estimate Std. Error z value        Pr(>|z|)    
ar1    1.120688   0.173833  6.4469 0.0000000001142 ***
ar2   -0.276165   0.152852 -1.8067          0.0708 .  
ma1   -0.054192   0.176020 -0.3079          0.7582    
drift  0.352255   0.354486  0.9937          0.3204    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
# drift is not significant and we ignore

Its information criteria is higher and some of coefficents are not significant. Still Arima(1,1,2) is best choice.

Checking models4 for y3

Code 
arima113 <- Arima(TMS.arimapart$y3,  # variable
                  order = c(1, 1, 3)  # (p,d,q) parameters
)
arima113
Series: TMS.arimapart$y3 
ARIMA(1,1,3) 

Coefficients:
         ar1     ma1     ma2      ma3
      0.8287  0.2506  0.0846  -0.1528
s.e.  0.0573  0.0845  0.0860   0.0800

sigma^2 = 1.031:  log likelihood = -427.66
AIC=865.32   AICc=865.53   BIC=883.82
Code 
coeftest(arima113)

z test of coefficients:

     Estimate Std. Error z value  Pr(>|z|)    
ar1  0.828733   0.057339 14.4533 < 2.2e-16 ***
ma1  0.250646   0.084506  2.9660  0.003017 ** 
ma2  0.084584   0.085995  0.9836  0.325316    
ma3 -0.152808   0.080038 -1.9092  0.056236 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
arima113_2 <- Arima(TMS.arimapart$y3,  # variable
                    order = c(1, 1, 3),  # (p,d,q) parameters
                    include.constant = TRUE)  # including a constant
arima113_2
Series: TMS.arimapart$y3 
ARIMA(1,1,3) with drift 

Coefficients:
         ar1    ma1     ma2      ma3   drift
      0.8198  0.257  0.0905  -0.1484  0.3524
s.e.  0.0597  0.086  0.0872   0.0806  0.3820

sigma^2 = 1.032:  log likelihood = -427.25
AIC=866.5   AICc=866.79   BIC=888.71
Code 
coeftest(arima113_2)

z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
ar1    0.819844   0.059687 13.7357 < 2.2e-16 ***
ma1    0.256971   0.086001  2.9880  0.002808 ** 
ma2    0.090468   0.087164  1.0379  0.299314    
ma3   -0.148404   0.080613 -1.8409  0.065631 .  
drift  0.352428   0.381977  0.9226  0.356194    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
# drift is not significant and we ignore it

Here Information criteria is lower than our previous last choice but some of coefficents are non-significant. For goodness-of-fit, I choose this model.

Checking models1 for y10

Code 
y10_arima212 <- Arima(TMS.arimapart$y3,  # variable
                  order = c(2, 1, 2)  # (p,d,q) parameters
)
y10_arima212
Series: TMS.arimapart$y3 
ARIMA(2,1,2) 

Coefficients:
         ar1     ar2     ma1     ma2
      0.3676  0.3115  0.7182  0.2992
s.e.  0.1843  0.1577  0.1789  0.0767

sigma^2 = 1.035:  log likelihood = -428.2
AIC=866.41   AICc=866.61   BIC=884.91
Code 
coeftest(y10_arima212)

z test of coefficients:

    Estimate Std. Error z value   Pr(>|z|)    
ar1 0.367610   0.184335  1.9942    0.04613 *  
ar2 0.311501   0.157742  1.9748    0.04830 *  
ma1 0.718177   0.178935  4.0136 0.00005979 ***
ma2 0.299187   0.076665  3.9025 0.00009520 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
y10_arima212_2 <- Arima(TMS.arimapart$y3,  # variable
                    order = c(2, 1, 2),  # (p,d,q) parameters
                    include.constant = TRUE)  # including a constant
y10_arima212_2
Series: TMS.arimapart$y3 
ARIMA(2,1,2) with drift 

Coefficients:
         ar1     ar2     ma1     ma2   drift
      0.3655  0.3040  0.7167  0.3015  0.3530
s.e.  0.1854  0.1578  0.1800  0.0769  0.3522

sigma^2 = 1.035:  log likelihood = -427.72
AIC=867.43   AICc=867.72   BIC=889.64
Code 
coeftest(y10_arima212_2)

z test of coefficients:

      Estimate Std. Error z value   Pr(>|z|)    
ar1   0.365550   0.185402  1.9717    0.04865 *  
ar2   0.304002   0.157799  1.9265    0.05404 .  
ma1   0.716699   0.179952  3.9827 0.00006813 ***
ma2   0.301517   0.076913  3.9203 0.00008846 ***
drift 0.352976   0.352158  1.0023    0.31619    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
 #We ignore drift

Our AIC information criteria is 866 and all coefficents are significant..

Checking models2 for y10

Code 
y10_arima313 <- Arima(TMS.arimapart$y3,  # variable
                      order = c(3, 1, 3)  # (p,d,q) parameters
)
y10_arima313
Series: TMS.arimapart$y3 
ARIMA(3,1,3) 

Coefficients:
          ar1     ar2     ar3     ma1     ma2     ma3
      -0.7283  0.3599  0.7221  1.8632  1.4738  0.3285
s.e.   0.0517  0.0603  0.0494  0.0691  0.0999  0.0667

sigma^2 = 1.016:  log likelihood = -425
AIC=863.99   AICc=864.38   BIC=889.9
Code 
coeftest(y10_arima313)

z test of coefficients:

     Estimate Std. Error  z value       Pr(>|z|)    
ar1 -0.728332   0.051680 -14.0930      < 2.2e-16 ***
ar2  0.359917   0.060277   5.9711 0.000000002357 ***
ar3  0.722057   0.049385  14.6209      < 2.2e-16 ***
ma1  1.863178   0.069140  26.9480      < 2.2e-16 ***
ma2  1.473794   0.099926  14.7489      < 2.2e-16 ***
ma3  0.328478   0.066700   4.9247 0.000000845005 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
summary(y10_arima313)
Series: TMS.arimapart$y3 
ARIMA(3,1,3) 

Coefficients:
          ar1     ar2     ar3     ma1     ma2     ma3
      -0.7283  0.3599  0.7221  1.8632  1.4738  0.3285
s.e.   0.0517  0.0603  0.0494  0.0691  0.0999  0.0667

sigma^2 = 1.016:  log likelihood = -425
AIC=863.99   AICc=864.38   BIC=889.9

Training set error measures:
                     ME      RMSE       MAE        MPE      MAPE      MASE
Training set 0.05203032 0.9961157 0.7770729 0.02783753 0.3367881 0.4403472
                   ACF1
Training set -0.0366561
Code 
#It is better than previous estimation since all coefficent are significant and It has lower information criteria.

It is better than previous estimation since all coefficents are significant and It has lower information criteria.

Checking models3 for y10

Code 
y10_arima314 <- Arima(TMS.arimapart$y3,  # variable
                      order = c(3, 1, 4)  # (p,d,q) parameters
)
y10_arima314
Series: TMS.arimapart$y3 
ARIMA(3,1,4) 

Coefficients:
          ar1     ar2     ar3     ma1     ma2     ma3      ma4
      -0.4776  0.6073  0.3532  1.5626  0.8918  0.1271  -0.0909
s.e.   0.2941  0.2483  0.2741  0.2973  0.4909  0.2589   0.1165

sigma^2 = 1.034:  log likelihood = -426.54
AIC=869.09   AICc=869.58   BIC=898.69
Code 
coeftest(y10_arima314)

z test of coefficients:

     Estimate Std. Error z value     Pr(>|z|)    
ar1 -0.477621   0.294107 -1.6240      0.10438    
ar2  0.607319   0.248328  2.4456      0.01446 *  
ar3  0.353157   0.274096  1.2884      0.19759    
ma1  1.562643   0.297296  5.2562 0.0000001471 ***
ma2  0.891780   0.490875  1.8167      0.06926 .  
ma3  0.127139   0.258862  0.4911      0.62332    
ma4 -0.090938   0.116510 -0.7805      0.43509    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
summary(y10_arima314)
Series: TMS.arimapart$y3 
ARIMA(3,1,4) 

Coefficients:
          ar1     ar2     ar3     ma1     ma2     ma3      ma4
      -0.4776  0.6073  0.3532  1.5626  0.8918  0.1271  -0.0909
s.e.   0.2941  0.2483  0.2741  0.2973  0.4909  0.2589   0.1165

sigma^2 = 1.034:  log likelihood = -426.54
AIC=869.09   AICc=869.58   BIC=898.69

Training set error measures:
                     ME     RMSE       MAE        MPE      MAPE      MASE
Training set 0.05503733 1.003119 0.7849284 0.02929993 0.3400635 0.4447987
                     ACF1
Training set -0.001214331

Arima(3,1,3) is better option because of all coefficents are significant and lower information criteria. After multiple times model estimation, my choice is ARIMA(1,1,3) for y3 and ARIMA(3,1,3) for y10.

et’s use this beautiful code taking into account of only information criteria.

Code 
arima.best.AIC <- 
  auto.arima(TMS.arimapart$y3,
             d = 1,             # parameter d of ARIMA model
             max.p = 6,         # Maximum value of p
             max.q = 6,         # Maximum value of q
             max.order = 12,    # maximum p+q
             start.p = 1,       # Starting value of p in stepwise procedure
             start.q = 1,       # Starting value of q in stepwise procedure
             ic = "aic",        # Information criterion to be used in model selection.
             stepwise = FALSE,  # if FALSE considers all models
             allowdrift = TRUE, # include a constant
             trace = TRUE)      # show summary of all models considered

 Fitting models using approximations to speed things up...

 ARIMA(0,1,0)                    : 1308.749
 ARIMA(0,1,0) with drift         : 1301.618
 ARIMA(0,1,1)                    : 1068.415
 ARIMA(0,1,1) with drift         : 1063.986
 ARIMA(0,1,2)                    : 933.5096
 ARIMA(0,1,2) with drift         : 931.3783
 ARIMA(0,1,3)                    : 901.523
 ARIMA(0,1,3) with drift         : 900.5292
 ARIMA(0,1,4)                    : 886.7096
 ARIMA(0,1,4) with drift         : 886.4316
 ARIMA(0,1,5)                    : 875.9129
 ARIMA(0,1,5) with drift         : 876.3092
 ARIMA(0,1,6)                    : 876.0923
 ARIMA(0,1,6) with drift         : 876.6186
 ARIMA(1,1,0)                    : 880.6256
 ARIMA(1,1,0) with drift         : 882.0434
 ARIMA(1,1,1)                    : 869.1987
 ARIMA(1,1,1) with drift         : 870.3328
 ARIMA(1,1,2)                    : 864.9552
 ARIMA(1,1,2) with drift         : 865.6984
 ARIMA(1,1,3)                    : 863.4044
 ARIMA(1,1,3) with drift         : 864.5231
 ARIMA(1,1,4)                    : 865.3663
 ARIMA(1,1,4) with drift         : 866.4425
 ARIMA(1,1,5)                    : 867.3638
 ARIMA(1,1,5) with drift         : 868.4157
 ARIMA(1,1,6)                    : 868.8967
 ARIMA(1,1,6) with drift         : 870.0535
 ARIMA(2,1,0)                    : 863.445
 ARIMA(2,1,0) with drift         : 864.6766
 ARIMA(2,1,1)                    : 865.3857
 ARIMA(2,1,1) with drift         : 866.6055
 ARIMA(2,1,2)                    : 861.7934
 ARIMA(2,1,2) with drift         : 862.7851
 ARIMA(2,1,3)                    : 861.6243
 ARIMA(2,1,3) with drift         : 862.9831
 ARIMA(2,1,4)                    : 863.3504
 ARIMA(2,1,4) with drift         : Inf
 ARIMA(2,1,5)                    : 865.2505
 ARIMA(2,1,5) with drift         : 866.9566
 ARIMA(2,1,6)                    : 866.5847
 ARIMA(2,1,6) with drift         : 867.9802
 ARIMA(3,1,0)                    : 866.2893
 ARIMA(3,1,0) with drift         : 867.4655
 ARIMA(3,1,1)                    : 866.2666
 ARIMA(3,1,1) with drift         : 867.3769
 ARIMA(3,1,2)                    : 864.7651
 ARIMA(3,1,2) with drift         : 865.7735
 ARIMA(3,1,3)                    : 863.9064
 ARIMA(3,1,3) with drift         : 865.4225
 ARIMA(3,1,4)                    : Inf
 ARIMA(3,1,4) with drift         : Inf
 ARIMA(3,1,5)                    : Inf
 ARIMA(3,1,5) with drift         : Inf
 ARIMA(3,1,6)                    : Inf
 ARIMA(3,1,6) with drift         : Inf
 ARIMA(4,1,0)                    : 862.4848
 ARIMA(4,1,0) with drift         : 863.9831
 ARIMA(4,1,1)                    : 864.269
 ARIMA(4,1,1) with drift         : 865.729
 ARIMA(4,1,2)                    : 865.7844
 ARIMA(4,1,2) with drift         : 867.1325
 ARIMA(4,1,3)                    : 865.5401
 ARIMA(4,1,3) with drift         : 867.0406
 ARIMA(4,1,4)                    : Inf
 ARIMA(4,1,4) with drift         : Inf
 ARIMA(4,1,5)                    : Inf
 ARIMA(4,1,5) with drift         : Inf
 ARIMA(4,1,6)                    : Inf
 ARIMA(4,1,6) with drift         : Inf
 ARIMA(5,1,0)                    : 864.1071
 ARIMA(5,1,0) with drift         : 865.569
 ARIMA(5,1,1)                    : 866.0995
 ARIMA(5,1,1) with drift         : 867.5647
 ARIMA(5,1,2)                    : 867.1588
 ARIMA(5,1,2) with drift         : 868.5158
 ARIMA(5,1,3)                    : Inf
 ARIMA(5,1,3) with drift         : Inf
 ARIMA(5,1,4)                    : 868.4973
 ARIMA(5,1,4) with drift         : 864.4925
 ARIMA(5,1,5)                    : 864.7485
 ARIMA(5,1,5) with drift         : 865.627
 ARIMA(5,1,6)                    : Inf
 ARIMA(5,1,6) with drift         : Inf
 ARIMA(6,1,0)                    : 866.8913
 ARIMA(6,1,0) with drift         : 868.3043
 ARIMA(6,1,1)                    : 866.0207
 ARIMA(6,1,1) with drift         : 867.4135
 ARIMA(6,1,2)                    : 867.6108
 ARIMA(6,1,2) with drift         : 868.9858
 ARIMA(6,1,3)                    : Inf
 ARIMA(6,1,3) with drift         : Inf
 ARIMA(6,1,4)                    : Inf
 ARIMA(6,1,4) with drift         : Inf
 ARIMA(6,1,5)                    : Inf
 ARIMA(6,1,5) with drift         : Inf
 ARIMA(6,1,6)                    : Inf
 ARIMA(6,1,6) with drift         : Inf

 Now re-fitting the best model(s) without approximations...




 Best model: ARIMA(2,1,3)
Code 
coeftest(arima.best.AIC)

z test of coefficients:

    Estimate Std. Error z value Pr(>|z|)
ar1  0.72485    0.45676  1.5869   0.1125
ar2  0.07962    0.34436  0.2312   0.8171
ma1  0.35309    0.45497  0.7761   0.4377
ma2  0.11944    0.18095  0.6601   0.5092
ma3 -0.13608    0.11440 -1.1895   0.2343

It has lowest information criteria but coefficents are not significant.

Code 
y10_arima.best.AIC <- auto.arima(TMS.arimapart$y10,
                                 d = 1,             # parameter d of ARIMA model
                                 max.p = 6,         # Maximum value of p
                                 max.q = 6,         # Maximum value of q
                                 max.order = 12,    # maximum p+q
                                 start.p = 1,       # Starting value of p in stepwise procedure
                                 start.q = 1,       # Starting value of q in stepwise procedure
                                 ic = "aic",        # Information criterion to be used in model selection.
                                 stepwise = FALSE,  # if FALSE considers all models
                                 allowdrift = TRUE, # include a constant
                                 trace = TRUE)      # show summary of all models considered

 Fitting models using approximations to speed things up...

 ARIMA(0,1,0)                    : 919.8482
 ARIMA(0,1,0) with drift         : 913.5013
 ARIMA(0,1,1)                    : 789.8636
 ARIMA(0,1,1) with drift         : 786.3076
 ARIMA(0,1,2)                    : 681.6209
 ARIMA(0,1,2) with drift         : 679.8767
 ARIMA(0,1,3)                    : 671.2678
 ARIMA(0,1,3) with drift         : 670.4239
 ARIMA(0,1,4)                    : 656.767
 ARIMA(0,1,4) with drift         : 656.5026
 ARIMA(0,1,5)                    : 650.1432
 ARIMA(0,1,5) with drift         : 650.2832
 ARIMA(0,1,6)                    : 650.3876
 ARIMA(0,1,6) with drift         : 650.6718
 ARIMA(1,1,0)                    : 659.6971
 ARIMA(1,1,0) with drift         : 660.6748
 ARIMA(1,1,1)                    : 653.0469
 ARIMA(1,1,1) with drift         : 654.3991
 ARIMA(1,1,2)                    : 647.6062
 ARIMA(1,1,2) with drift         : 648.5343
 ARIMA(1,1,3)                    : 640.7574
 ARIMA(1,1,3) with drift         : 642.1934
 ARIMA(1,1,4)                    : 640.6352
 ARIMA(1,1,4) with drift         : 641.8894
 ARIMA(1,1,5)                    : 642.4356
 ARIMA(1,1,5) with drift         : 643.7541
 ARIMA(1,1,6)                    : 644.2745
 ARIMA(1,1,6) with drift         : 645.6286
 ARIMA(2,1,0)                    : 648.6368
 ARIMA(2,1,0) with drift         : 650.1194
 ARIMA(2,1,1)                    : 643.3463
 ARIMA(2,1,1) with drift         : 644.7078
 ARIMA(2,1,2)                    : 642.3904
 ARIMA(2,1,2) with drift         : 643.5202
 ARIMA(2,1,3)                    : 640.711
 ARIMA(2,1,3) with drift         : 642.1819
 ARIMA(2,1,4)                    : 641.8232
 ARIMA(2,1,4) with drift         : 643.3761
 ARIMA(2,1,5)                    : 637.2944
 ARIMA(2,1,5) with drift         : 638.3904
 ARIMA(2,1,6)                    : 636.8118
 ARIMA(2,1,6) with drift         : 638.2124
 ARIMA(3,1,0)                    : 643.2229
 ARIMA(3,1,0) with drift         : 644.5093
 ARIMA(3,1,1)                    : 644.204
 ARIMA(3,1,1) with drift         : 645.4756
 ARIMA(3,1,2)                    : 646.1067
 ARIMA(3,1,2) with drift         : 647.1282
 ARIMA(3,1,3)                    : 643.1918
 ARIMA(3,1,3) with drift         : 644.7982
 ARIMA(3,1,4)                    : 643.5063
 ARIMA(3,1,4) with drift         : Inf
 ARIMA(3,1,5)                    : 646.5962
 ARIMA(3,1,5) with drift         : 647.5236
 ARIMA(3,1,6)                    : 647.2481
 ARIMA(3,1,6) with drift         : 648.4942
 ARIMA(4,1,0)                    : 639.3739
 ARIMA(4,1,0) with drift         : 640.9589
 ARIMA(4,1,1)                    : 641.2724
 ARIMA(4,1,1) with drift         : 642.8819
 ARIMA(4,1,2)                    : 641.4201
 ARIMA(4,1,2) with drift         : 642.8553
 ARIMA(4,1,3)                    : 639.8416
 ARIMA(4,1,3) with drift         : 641.4752
 ARIMA(4,1,4)                    : 641.8407
 ARIMA(4,1,4) with drift         : 643.4752
 ARIMA(4,1,5)                    : 643.7579
 ARIMA(4,1,5) with drift         : 645.3771
 ARIMA(4,1,6)                    : 640.2968
 ARIMA(4,1,6) with drift         : 641.9116
 ARIMA(5,1,0)                    : 641.6887
 ARIMA(5,1,0) with drift         : 643.3181
 ARIMA(5,1,1)                    : Inf
 ARIMA(5,1,1) with drift         : Inf
 ARIMA(5,1,2)                    : 641.875
 ARIMA(5,1,2) with drift         : 643.3801
 ARIMA(5,1,3)                    : 642.7136
 ARIMA(5,1,3) with drift         : 644.3284
 ARIMA(5,1,4)                    : 638.5627
 ARIMA(5,1,4) with drift         : 639.7602
 ARIMA(5,1,5)                    : 640.5477
 ARIMA(5,1,5) with drift         : 641.7155
 ARIMA(5,1,6)                    : 642.282
 ARIMA(5,1,6) with drift         : 643.562
 ARIMA(6,1,0)                    : 642.3694
 ARIMA(6,1,0) with drift         : 643.8719
 ARIMA(6,1,1)                    : 642.9665
 ARIMA(6,1,1) with drift         : 644.5178
 ARIMA(6,1,2)                    : 642.1514
 ARIMA(6,1,2) with drift         : 643.9489
 ARIMA(6,1,3)                    : 637.0177
 ARIMA(6,1,3) with drift         : 638.7387
 ARIMA(6,1,4)                    : Inf
 ARIMA(6,1,4) with drift         : Inf
 ARIMA(6,1,5)                    : Inf
 ARIMA(6,1,5) with drift         : Inf
 ARIMA(6,1,6)                    : Inf
 ARIMA(6,1,6) with drift         : Inf

 Now re-fitting the best model(s) without approximations...

 ARIMA(0,1,0)                    : 920.0809
 ARIMA(0,1,0) with drift         : 913.734
 ARIMA(0,1,1)                    : 790.312
 ARIMA(0,1,1) with drift         : 786.7002
 ARIMA(0,1,2)                    : 682.6714
 ARIMA(0,1,2) with drift         : 680.9753
 ARIMA(0,1,3)                    : 672.4007
 ARIMA(0,1,3) with drift         : 671.4942
 ARIMA(0,1,4)                    : 657.8685
 ARIMA(0,1,4) with drift         : 657.5873
 ARIMA(0,1,5)                    : 651.2064
 ARIMA(0,1,5) with drift         : 651.3528
 ARIMA(0,1,6)                    : 651.4598
 ARIMA(0,1,6) with drift         : 651.7762
 ARIMA(1,1,0)                    : 659.8061
 ARIMA(1,1,0) with drift         : 660.7804
 ARIMA(1,1,1)                    : 653.162
 ARIMA(1,1,1) with drift         : 654.4592
 ARIMA(1,1,2)                    : 647.7939
 ARIMA(1,1,2) with drift         : 648.7721
 ARIMA(1,1,3)                    : 640.9856
 ARIMA(1,1,3) with drift         : 642.2947
 ARIMA(1,1,4)                    : 640.8259
 ARIMA(1,1,4) with drift         : 642.0333
 ARIMA(1,1,5)                    : 642.62
 ARIMA(1,1,5) with drift         : 643.865
 ARIMA(1,1,6)                    : 644.4629
 ARIMA(1,1,6) with drift         : 645.7222
 ARIMA(2,1,0)                    : 649.6427
 ARIMA(2,1,0) with drift         : 650.9792
 ARIMA(2,1,1)                    : 644.2041
 ARIMA(2,1,1) with drift         : 645.4783
 ARIMA(2,1,2)                    : 642.5037
 ARIMA(2,1,2) with drift         : 643.5677
 ARIMA(2,1,3)                    : 641.0608
 ARIMA(2,1,3) with drift         : 642.3346
 ARIMA(2,1,4)                    : 642.7189
 ARIMA(2,1,4) with drift         : 643.943
 ARIMA(2,1,5)                    : Inf
 ARIMA(2,1,5) with drift         : Inf
 ARIMA(2,1,6)                    : Inf
 ARIMA(2,1,6) with drift         : Inf
 ARIMA(3,1,0)                    : 643.784
 ARIMA(3,1,0) with drift         : 644.872
 ARIMA(3,1,1)                    : 644.1032
 ARIMA(3,1,1) with drift         : 645.2326
 ARIMA(3,1,2)                    : 637.2385
 ARIMA(3,1,2) with drift         : 638.365
 ARIMA(3,1,3)                    : 638.2398
 ARIMA(3,1,3) with drift         : 639.4768
 ARIMA(3,1,4)                    : 643.055
 ARIMA(3,1,4) with drift         : 640.9909
 ARIMA(3,1,5)                    : 639.926
 ARIMA(3,1,5) with drift         : 641.224
 ARIMA(3,1,6)                    : 641.8064
 ARIMA(3,1,6) with drift         : 643.0855
 ARIMA(4,1,0)                    : 642.328
 ARIMA(4,1,0) with drift         : 643.5732
 ARIMA(4,1,1)                    : 644.0604
 ARIMA(4,1,1) with drift         : 645.3405
 ARIMA(4,1,2)                    : 637.9736
 ARIMA(4,1,2) with drift         : 639.2215
 ARIMA(4,1,3)                    : Inf
 ARIMA(4,1,3) with drift         : Inf
 ARIMA(4,1,4)                    : 642.2376
 ARIMA(4,1,4) with drift         : Inf
 ARIMA(4,1,5)                    : 641.8622
 ARIMA(4,1,5) with drift         : 643.1605
 ARIMA(4,1,6)                    : Inf
 ARIMA(4,1,6) with drift         : Inf
 ARIMA(5,1,0)                    : 643.6852
 ARIMA(5,1,0) with drift         : 644.9918
 ARIMA(5,1,1)                    : Inf
 ARIMA(5,1,1) with drift         : Inf
 ARIMA(5,1,2)                    : 639.3584
 ARIMA(5,1,2) with drift         : 640.5435
 ARIMA(5,1,3)                    : 641.3514
 ARIMA(5,1,3) with drift         : 642.4405
 ARIMA(5,1,4)                    : Inf
 ARIMA(5,1,4) with drift         : Inf
 ARIMA(5,1,5)                    : Inf
 ARIMA(5,1,5) with drift         : Inf
 ARIMA(5,1,6)                    : Inf
 ARIMA(5,1,6) with drift         : Inf
 ARIMA(6,1,0)                    : 643.4006
 ARIMA(6,1,0) with drift         : 644.5821
 ARIMA(6,1,1)                    : 643.9905
 ARIMA(6,1,1) with drift         : 645.205
 ARIMA(6,1,2)                    : 640.9941
 ARIMA(6,1,2) with drift         : 642.2265
 ARIMA(6,1,3)                    : Inf
 ARIMA(6,1,3) with drift         : Inf
 ARIMA(6,1,4)                    : Inf
 ARIMA(6,1,4) with drift         : Inf
 ARIMA(6,1,5)                    : Inf
 ARIMA(6,1,5) with drift         : Inf
 ARIMA(6,1,6)                    : Inf
 ARIMA(6,1,6) with drift         : Inf





 Best model: ARIMA(3,1,2)
Code 
coeftest(y10_arima.best.AIC)

z test of coefficients:

     Estimate Std. Error z value           Pr(>|z|)    
ar1 -0.449081   0.105266 -4.2662 0.0000198867375112 ***
ar2  0.303120   0.100606  3.0129           0.002587 ** 
ar3  0.536404   0.089012  6.0262 0.0000000016789571 ***
ma1  1.144355   0.089870 12.7334          < 2.2e-16 ***
ma2  0.757863   0.101857  7.4405 0.0000000000001003 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

It has lowest information criteria and all of the coefficents are significant. For y10, I will use ARIMA(3,1,2) and for y3, I will use ARIMA(1,1,2)

Model Diagnostics for y3

Plotting residuals

Code 
y3_arima <- arima113
plot(resid(y3_arima))

Code 
tibble(
  date = index(TMS.arimapart),
  resid = arima113 %>% resid() %>% as.numeric()
) %>%
  ggplot(aes(date, resid)) +
  geom_line(col = "royalblue3") +
  theme_bw()

Visually we can conclude that our residuals are whit noise.

ACF and PACF plot of residuals

Code 
par(mfrow = c(2, 1)) 
acf(resid(y3_arima), 
    lag.max = 36,
    ylim = c(-0.5, 0.5), 
    lwd = 5, col = "dark green",
    na.action = na.pass)
pacf(resid(y3_arima), 
     lag.max = 36, 
     lwd = 5, col = "dark green",
     na.action = na.pass)

Code 
par(mfrow = c(1, 1))

Visually till high lags, we cannot see autocorrelation, however it should be better to test statistically.

Ljung box test for residuals

Code 
Box.test(resid(y3_arima), type = "Ljung-Box", lag = 5)

    Box-Ljung test

data:  resid(y3_arima)
X-squared = 0.028645, df = 5, p-value = 1
Code 
Box.test(resid(y3_arima), type = "Ljung-Box", lag = 10)

    Box-Ljung test

data:  resid(y3_arima)
X-squared = 3.0289, df = 10, p-value = 0.9807
Code 
Box.test(resid(y3_arima), type = "Ljung-Box", lag = 15)

    Box-Ljung test

data:  resid(y3_arima)
X-squared = 5.5429, df = 15, p-value = 0.9864

P values are higher than significance level, so we cannot reject our residuals are not autocorrelated. Residuals are free from autocorrelation till lag=15, we can use it for forecasting.

Model Diagnostics for y10(plotting residuals)

Code 
y10_arima <- y10_arima.best.AIC
plot(resid(y10_arima))

Code 
tibble(
  date = index(TMS.arimapart),
  resid = y10_arima %>% resid() %>% as.numeric()
) %>%
  ggplot(aes(date, resid)) +
  geom_line(col = "royalblue3") +
  theme_bw()

ACF and PACF plot of residuals

Code 
par(mfrow = c(2, 1)) 
acf(resid(y10_arima), 
    lag.max = 36,
    ylim = c(-0.5, 0.5), 
    lwd = 5, col = "dark green",
    na.action = na.pass)
pacf(resid(y10_arima), 
     lag.max = 36, 
     lwd = 5, col = "dark green",
     na.action = na.pass)

Code 
par(mfrow = c(1, 1))

Our residuals are white noise but still we need to test statistically.

Ljung box test for residuals

Code 
Box.test(resid(y10_arima), type = "Ljung-Box", lag = 5)

    Box-Ljung test

data:  resid(y10_arima)
X-squared = 0.21602, df = 5, p-value = 0.9989
Code 
Box.test(resid(y10_arima), type = "Ljung-Box", lag = 10)

    Box-Ljung test

data:  resid(y10_arima)
X-squared = 1.1619, df = 10, p-value = 0.9997
Code 
Box.test(resid(y10_arima), type = "Ljung-Box", lag = 15)

    Box-Ljung test

data:  resid(y10_arima)
X-squared = 6.363, df = 15, p-value = 0.973

Residuals of y10_arima are free from autocorrelation

Forecasting for y3

Code 
tail(TMS.arimapart, 20)
               y3    y10   dy3  dy10
2022-03-08 302.89 253.33 -3.87 -1.81
2022-03-09 300.11 251.60 -2.78 -1.73
2022-03-10 298.22 250.84 -1.89 -0.76
2022-03-11 297.11 250.14 -1.11 -0.70
2022-03-12 296.66 249.92 -0.45 -0.22
2022-03-13 296.83 250.31  0.17  0.39
2022-03-14 296.04 249.75 -0.79 -0.56
2022-03-15 293.82 249.04 -2.22 -0.71
2022-03-16 291.49 247.09 -2.33 -1.95
2022-03-17 290.05 246.74 -1.44 -0.35
2022-03-18 288.80 246.23 -1.25 -0.51
2022-03-19 288.38 246.42 -0.42  0.19
2022-03-20 288.44 246.38  0.06 -0.04
2022-03-21 288.86 246.68  0.42  0.30
2022-03-22 289.20 246.30  0.34 -0.38
2022-03-23 288.83 245.85 -0.37 -0.45
2022-03-24 287.98 245.64 -0.85 -0.21
2022-03-25 287.54 245.77 -0.44  0.13
2022-03-26 286.01 244.87 -1.53 -0.90
2022-03-27 285.43 244.08 -0.58 -0.79
Code 
forecasts <- forecast(y3_arima, # model for prediction
                      h = 20) # how many periods outside the sample
forecasts
    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
301       285.0073 283.7059 286.3087 283.0169 286.9976
302       284.9394 281.9366 287.9422 280.3470 289.5318
303       284.7407 279.7547 289.7267 277.1153 292.3662
304       284.5761 277.6291 291.5232 273.9515 295.2007
305       284.4397 275.5493 293.3301 270.8429 298.0365
306       284.3266 273.5240 295.1293 267.8055 300.8478
307       284.2329 271.5610 296.9049 264.8529 303.6130
308       284.1553 269.6646 298.6460 261.9937 306.3169
309       284.0909 267.8366 300.3453 259.2320 308.9498
310       284.0376 266.0767 301.9985 256.5688 311.5065
311       283.9934 264.3833 303.6035 254.0024 313.9844
312       283.9568 262.7541 305.1595 251.5301 316.3835
313       283.9264 261.1860 306.6668 249.1480 318.7049
314       283.9013 259.6759 308.1267 246.8518 320.9508
315       283.8804 258.2204 309.5405 244.6368 323.1240
316       283.8632 256.8163 310.9100 242.4985 325.2278
317       283.8488 255.4603 312.2374 240.4323 327.2654
318       283.8370 254.1493 313.5247 238.4336 329.2404
319       283.8271 252.8804 314.7738 236.4983 331.1560
320       283.8190 251.6510 315.9870 234.6223 333.0157
Code 
str(forecasts)
List of 10
 $ method   : chr "ARIMA(1,1,3)"
 $ model    :List of 18
  ..$ coef     : Named num [1:4] 0.8287 0.2506 0.0846 -0.1528
  .. ..- attr(*, "names")= chr [1:4] "ar1" "ma1" "ma2" "ma3"
  ..$ sigma2   : num 1.03
  ..$ var.coef : num [1:4, 1:4] 0.00329 -0.00355 -0.00364 -0.00302 -0.00355 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:4] "ar1" "ma1" "ma2" "ma3"
  .. .. ..$ : chr [1:4] "ar1" "ma1" "ma2" "ma3"
  ..$ mask     : logi [1:4] TRUE TRUE TRUE TRUE
  ..$ loglik   : num -428
  ..$ aic      : num 865
  ..$ arma     : int [1:7] 1 3 0 0 1 1 0
  ..$ residuals: Time-Series [1:300] from 1 to 300: 0.17357 -0.000321 2.04479 -0.088908 1.294769 ...
  ..$ call     : language Arima(y = TMS.arimapart$y3, order = c(1, 1, 3))
  ..$ series   : chr "TMS.arimapart$y3"
  ..$ code     : int 0
  ..$ n.cond   : int 0
  ..$ nobs     : int 299
  ..$ model    :List of 10
  .. ..$ phi  : num 0.829
  .. ..$ theta: num [1:3] 0.2506 0.0846 -0.1528
  .. ..$ Delta: num 1
  .. ..$ Z    : num [1:5] 1 0 0 0 1
  .. ..$ a    : num [1:5] -0.58 0.058 0.282 -0.142 286.01
  .. ..$ P    : num [1:5, 1:5] 0.00 0.00 0.00 0.00 -1.52e-17 ...
  .. ..$ T    : num [1:5, 1:5] 0.829 0 0 0 1 ...
  .. ..$ V    : num [1:5, 1:5] 1 0.2506 0.0846 -0.1528 0 ...
  .. ..$ h    : num 0
  .. ..$ Pn   : num [1:5, 1:5] 1.00 2.51e-01 8.46e-02 -1.53e-01 1.30e-17 ...
  ..$ aicc     : num 866
  ..$ bic      : num 884
  ..$ x        : Time-Series [1:300] from 1 to 300: 174 174 176 178 181 ...
  ..$ fitted   : Time-Series [1:300] from 1 to 300: 173 174 174 178 180 ...
  ..- attr(*, "class")= chr [1:3] "forecast_ARIMA" "ARIMA" "Arima"
 $ level    : num [1:2] 80 95
 $ mean     : Time-Series [1:20] from 301 to 320: 285 285 285 285 284 ...
 $ lower    : Time-Series [1:20, 1:2] from 301 to 320: 284 282 280 278 276 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "80%" "95%"
 $ upper    : Time-Series [1:20, 1:2] from 301 to 320: 286 288 290 292 293 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "80%" "95%"
 $ x        : Time-Series [1:300] from 1 to 300: 174 174 176 178 181 ...
 $ series   : chr "TMS.arimapart$y3"
 $ fitted   : Time-Series [1:300] from 1 to 300: 173 174 174 178 180 ...
 $ residuals: Time-Series [1:300] from 1 to 300: 0.17357 -0.000321 2.04479 -0.088908 1.294769 ...
 - attr(*, "class")= chr "forecast"
Code 
forecasts$mean
Time Series:
Start = 301 
End = 320 
Frequency = 1 
 [1] 285.0073 284.9394 284.7407 284.5761 284.4397 284.3266 284.2329 284.1553
 [9] 284.0909 284.0376 283.9934 283.9568 283.9264 283.9013 283.8804 283.8632
[17] 283.8488 283.8370 283.8271 283.8190
Code 
class(forecasts$mean)
[1] "ts"
Code 
as.numeric(forecasts$mean)
 [1] 285.0073 284.9394 284.7407 284.5761 284.4397 284.3266 284.2329 284.1553
 [9] 284.0909 284.0376 283.9934 283.9568 283.9264 283.9013 283.8804 283.8632
[17] 283.8488 283.8370 283.8271 283.8190
Code 
forecasts$lower
Time Series:
Start = 301 
End = 320 
Frequency = 1 
         80%      95%
301 283.7059 283.0169
302 281.9366 280.3470
303 279.7547 277.1153
304 277.6291 273.9515
305 275.5493 270.8429
306 273.5240 267.8055
307 271.5610 264.8529
308 269.6646 261.9937
309 267.8366 259.2320
310 266.0767 256.5688
311 264.3833 254.0024
312 262.7541 251.5301
313 261.1860 249.1480
314 259.6759 246.8518
315 258.2204 244.6368
316 256.8163 242.4985
317 255.4603 240.4323
318 254.1493 238.4336
319 252.8804 236.4983
320 251.6510 234.6223
Code 
forecasts$upper
Time Series:
Start = 301 
End = 320 
Frequency = 1 
         80%      95%
301 286.3087 286.9976
302 287.9422 289.5318
303 289.7267 292.3662
304 291.5232 295.2007
305 293.3301 298.0365
306 295.1293 300.8478
307 296.9049 303.6130
308 298.6460 306.3169
309 300.3453 308.9498
310 301.9985 311.5065
311 303.6035 313.9844
312 305.1595 316.3835
313 306.6668 318.7049
314 308.1267 320.9508
315 309.5405 323.1240
316 310.9100 325.2278
317 312.2374 327.2654
318 313.5247 329.2404
319 314.7738 331.1560
320 315.9870 333.0157
Code 
forecasts_data <- data.frame(f_mean  = as.numeric(forecasts$mean),
                             f_lower = as.numeric(forecasts$lower[, 2]),
                             f_upper = as.numeric(forecasts$upper[, 2]))
forecasts_data
f_mean f_lower f_upper
285.0073 283.0169 286.9976
284.9394 280.3470 289.5318
284.7407 277.1153 292.3662
284.5761 273.9515 295.2007
284.4397 270.8429 298.0365
284.3266 267.8055 300.8478
284.2329 264.8529 303.6130
284.1553 261.9937 306.3169
284.0909 259.2320 308.9498
284.0376 256.5688 311.5065
283.9934 254.0024 313.9844
283.9568 251.5301 316.3835
283.9264 249.1480 318.7049
283.9013 246.8518 320.9508
283.8804 244.6368 323.1240
283.8632 242.4985 325.2278
283.8488 240.4323 327.2654
283.8370 238.4336 329.2404
283.8271 236.4983 331.1560
283.8190 234.6223 333.0157
Code 
forecasts_data_xts_y3 <- TMS.arimapart[281:300, "y3"]
forecasts_xts_1 <- xts(forecasts_data,
                       order.by = index(forecasts_data_xts_y3))
forecasts_xts_1
             f_mean  f_lower  f_upper
2022-03-08 285.0073 283.0169 286.9976
2022-03-09 284.9394 280.3470 289.5318
2022-03-10 284.7407 277.1153 292.3662
2022-03-11 284.5761 273.9515 295.2007
2022-03-12 284.4397 270.8429 298.0365
2022-03-13 284.3266 267.8055 300.8478
2022-03-14 284.2329 264.8529 303.6130
2022-03-15 284.1553 261.9937 306.3169
2022-03-16 284.0909 259.2320 308.9498
2022-03-17 284.0376 256.5688 311.5065
2022-03-18 283.9934 254.0024 313.9844
2022-03-19 283.9568 251.5301 316.3835
2022-03-20 283.9264 249.1480 318.7049
2022-03-21 283.9013 246.8518 320.9508
2022-03-22 283.8804 244.6368 323.1240
2022-03-23 283.8632 242.4985 325.2278
2022-03-24 283.8488 240.4323 327.2654
2022-03-25 283.8370 238.4336 329.2404
2022-03-26 283.8271 236.4983 331.1560
2022-03-27 283.8190 234.6223 333.0157
Code 
forecasting_y3 <- merge(financial_data$y3, forecasts_xts_1)
tail(forecasting_y3, n = 20)
               y3   f_mean  f_lower  f_upper
2022-03-08 302.89 285.0073 283.0169 286.9976
2022-03-09 300.11 284.9394 280.3470 289.5318
2022-03-10 298.22 284.7407 277.1153 292.3662
2022-03-11 297.11 284.5761 273.9515 295.2007
2022-03-12 296.66 284.4397 270.8429 298.0365
2022-03-13 296.83 284.3266 267.8055 300.8478
2022-03-14 296.04 284.2329 264.8529 303.6130
2022-03-15 293.82 284.1553 261.9937 306.3169
2022-03-16 291.49 284.0909 259.2320 308.9498
2022-03-17 290.05 284.0376 256.5688 311.5065
2022-03-18 288.80 283.9934 254.0024 313.9844
2022-03-19 288.38 283.9568 251.5301 316.3835
2022-03-20 288.44 283.9264 249.1480 318.7049
2022-03-21 288.86 283.9013 246.8518 320.9508
2022-03-22 289.20 283.8804 244.6368 323.1240
2022-03-23 288.83 283.8632 242.4985 325.2278
2022-03-24 287.98 283.8488 240.4323 327.2654
2022-03-25 287.54 283.8370 238.4336 329.2404
2022-03-26 286.01 283.8271 236.4983 331.1560
2022-03-27 285.43 283.8190 234.6223 333.0157

We only need 95 level confidence iterval, so we removed 80 upper and lower boundaries.

Plotting forecast of y3

Code 
par(mfrow = c(1, 1))

plot(forecasting_y3[, c("y3", "f_mean", "f_lower", "f_upper")], 
     major.ticks = "years", 
     grid.ticks.on = "years",
     grid.ticks.lty = 3,
     main = "20 day forecast of y3",
     col = c("black", "blue", "red", "red"))

We need to dive deeper

Code 
plot(forecasting_y3["2022-02/", # limit the rows
           c("y3", "f_mean", "f_lower", "f_upper")], 
     major.ticks = "days", 
     grid.ticks.on = "days",
     grid.ticks.lty = 3,
     main = "20 day forecast of y3",
     col = c("black", "blue", "red", "red"))

We can see that our blue line(our forecasts) is located between red lines(boundaries). In the first days of comparisons, deviation from original series is high. However, it approachs black line in the end of sample.

Forecasting for y10

Code 
tail(TMS.arimapart, 20)
               y3    y10   dy3  dy10
2022-03-08 302.89 253.33 -3.87 -1.81
2022-03-09 300.11 251.60 -2.78 -1.73
2022-03-10 298.22 250.84 -1.89 -0.76
2022-03-11 297.11 250.14 -1.11 -0.70
2022-03-12 296.66 249.92 -0.45 -0.22
2022-03-13 296.83 250.31  0.17  0.39
2022-03-14 296.04 249.75 -0.79 -0.56
2022-03-15 293.82 249.04 -2.22 -0.71
2022-03-16 291.49 247.09 -2.33 -1.95
2022-03-17 290.05 246.74 -1.44 -0.35
2022-03-18 288.80 246.23 -1.25 -0.51
2022-03-19 288.38 246.42 -0.42  0.19
2022-03-20 288.44 246.38  0.06 -0.04
2022-03-21 288.86 246.68  0.42  0.30
2022-03-22 289.20 246.30  0.34 -0.38
2022-03-23 288.83 245.85 -0.37 -0.45
2022-03-24 287.98 245.64 -0.85 -0.21
2022-03-25 287.54 245.77 -0.44  0.13
2022-03-26 286.01 244.87 -1.53 -0.90
2022-03-27 285.43 244.08 -0.58 -0.79
Code 
forecasts <- forecast(y10_arima, # model for prediction
                      h = 20) # how many periods outside the sample
forecasts
    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
301       243.2799 242.3921 244.1677 241.9221 244.6377
302       242.7846 241.0371 244.5320 240.1121 245.4571
303       242.3408 239.5547 245.1268 238.0799 246.6016
304       241.9608 238.1943 245.7272 236.2004 247.7211
305       241.7312 236.9755 246.4868 234.4580 249.0043
306       241.4810 235.7462 247.2159 232.7103 250.2517
307       241.3199 234.6647 247.9752 231.1416 251.4983
308       241.1933 233.6320 248.7546 229.6292 252.7574
309       241.0672 232.6457 249.4886 228.1877 253.9467
310       240.9990 231.7554 250.2426 226.8621 255.1359
311       240.9235 230.8852 250.9617 225.5713 256.2756
312       240.8691 230.0770 251.6612 224.3640 257.3742
313       240.8340 229.3148 252.3532 223.2169 258.4511
314       240.7928 228.5761 253.0094 222.1090 259.4765
315       240.7715 227.8865 253.6565 221.0655 260.4775
316       240.7497 227.2187 254.2808 220.0558 261.4437
317       240.7309 226.5791 254.8828 219.0875 262.3743
318       240.7214 225.9696 255.4732 218.1604 263.2823
319       240.7083 225.3760 256.0406 217.2595 264.1571
320       240.7012 224.8077 256.5946 216.3942 265.0081
Code 
str(forecasts)
List of 10
 $ method   : chr "ARIMA(3,1,2)"
 $ model    :List of 18
  ..$ coef     : Named num [1:5] -0.449 0.303 0.536 1.144 0.758
  .. ..- attr(*, "names")= chr [1:5] "ar1" "ar2" "ar3" "ma1" ...
  ..$ sigma2   : num 0.48
  ..$ var.coef : num [1:5, 1:5] 0.01108 -0.00075 -0.00582 -0.00803 -0.00281 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:5] "ar1" "ar2" "ar3" "ma1" ...
  .. .. ..$ : chr [1:5] "ar1" "ar2" "ar3" "ma1" ...
  ..$ mask     : logi [1:5] TRUE TRUE TRUE TRUE TRUE
  ..$ loglik   : num -313
  ..$ aic      : num 637
  ..$ arma     : int [1:7] 3 2 0 0 1 1 0
  ..$ residuals: Time-Series [1:300] from 1 to 300: 0.188 0.012 0.926 0.456 1.237 ...
  ..$ call     : language auto.arima(y = TMS.arimapart$y10, d = 1, max.p = 6, max.q = 6, max.order = 12,      start.p = 1, start.q = 1, ic | __truncated__ ...
  ..$ series   : chr "TMS.arimapart$y10"
  ..$ code     : int 0
  ..$ n.cond   : int 0
  ..$ nobs     : int 299
  ..$ model    :List of 10
  .. ..$ phi  : num [1:3] -0.449 0.303 0.536
  .. ..$ theta: num [1:2] 1.144 0.758
  .. ..$ Delta: num 1
  .. ..$ Z    : num [1:4] 1 0 0 1
  .. ..$ a    : num [1:4] -0.79 -1.155 -0.615 244.87
  .. ..$ P    : num [1:4, 1:4] -1.11e-16 0.00 0.00 4.99e-17 0.00 ...
  .. ..$ T    : num [1:4, 1:4] -0.449 0.303 0.536 1 1 ...
  .. ..$ V    : num [1:4, 1:4] 1 1.144 0.758 0 1.144 ...
  .. ..$ h    : num 0
  .. ..$ Pn   : num [1:4, 1:4] 1.00 1.14 7.58e-01 7.83e-17 1.14 ...
  ..$ bic      : num 659
  ..$ aicc     : num 638
  ..$ x        : Time-Series [1:300] from 1 to 300: 188 188 189 190 193 ...
  ..$ fitted   : Time-Series [1:300] from 1 to 300: 188 188 188 190 191 ...
  ..- attr(*, "class")= chr [1:3] "forecast_ARIMA" "ARIMA" "Arima"
 $ level    : num [1:2] 80 95
 $ mean     : Time-Series [1:20] from 301 to 320: 243 243 242 242 242 ...
 $ lower    : Time-Series [1:20, 1:2] from 301 to 320: 242 241 240 238 237 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "80%" "95%"
 $ upper    : Time-Series [1:20, 1:2] from 301 to 320: 244 245 245 246 246 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "80%" "95%"
 $ x        : Time-Series [1:300] from 1 to 300: 188 188 189 190 193 ...
 $ series   : chr "TMS.arimapart$y10"
 $ fitted   : Time-Series [1:300] from 1 to 300: 188 188 188 190 191 ...
 $ residuals: Time-Series [1:300] from 1 to 300: 0.188 0.012 0.926 0.456 1.237 ...
 - attr(*, "class")= chr "forecast"
Code 
forecasts$mean
Time Series:
Start = 301 
End = 320 
Frequency = 1 
 [1] 243.2799 242.7846 242.3408 241.9608 241.7312 241.4810 241.3199 241.1933
 [9] 241.0672 240.9990 240.9235 240.8691 240.8340 240.7928 240.7715 240.7497
[17] 240.7309 240.7214 240.7083 240.7012
Code 
class(forecasts$mean)
[1] "ts"
Code 
as.numeric(forecasts$mean)
 [1] 243.2799 242.7846 242.3408 241.9608 241.7312 241.4810 241.3199 241.1933
 [9] 241.0672 240.9990 240.9235 240.8691 240.8340 240.7928 240.7715 240.7497
[17] 240.7309 240.7214 240.7083 240.7012
Code 
forecasts$lower
Time Series:
Start = 301 
End = 320 
Frequency = 1 
         80%      95%
301 242.3921 241.9221
302 241.0371 240.1121
303 239.5547 238.0799
304 238.1943 236.2004
305 236.9755 234.4580
306 235.7462 232.7103
307 234.6647 231.1416
308 233.6320 229.6292
309 232.6457 228.1877
310 231.7554 226.8621
311 230.8852 225.5713
312 230.0770 224.3640
313 229.3148 223.2169
314 228.5761 222.1090
315 227.8865 221.0655
316 227.2187 220.0558
317 226.5791 219.0875
318 225.9696 218.1604
319 225.3760 217.2595
320 224.8077 216.3942
Code 
forecasts$upper
Time Series:
Start = 301 
End = 320 
Frequency = 1 
         80%      95%
301 244.1677 244.6377
302 244.5320 245.4571
303 245.1268 246.6016
304 245.7272 247.7211
305 246.4868 249.0043
306 247.2159 250.2517
307 247.9752 251.4983
308 248.7546 252.7574
309 249.4886 253.9467
310 250.2426 255.1359
311 250.9617 256.2756
312 251.6612 257.3742
313 252.3532 258.4511
314 253.0094 259.4765
315 253.6565 260.4775
316 254.2808 261.4437
317 254.8828 262.3743
318 255.4732 263.2823
319 256.0406 264.1571
320 256.5946 265.0081
Code 
forecasts_data <- data.frame(f_mean  = as.numeric(forecasts$mean),
                             f_lower = as.numeric(forecasts$lower[, 2]),
                             f_upper = as.numeric(forecasts$upper[, 2]))
forecasts_data
f_mean f_lower f_upper
243.2799 241.9221 244.6377
242.7846 240.1121 245.4571
242.3408 238.0799 246.6016
241.9608 236.2004 247.7211
241.7312 234.4580 249.0043
241.4810 232.7103 250.2517
241.3199 231.1416 251.4983
241.1933 229.6292 252.7574
241.0672 228.1877 253.9467
240.9990 226.8621 255.1359
240.9235 225.5713 256.2756
240.8691 224.3640 257.3742
240.8340 223.2169 258.4511
240.7928 222.1090 259.4765
240.7715 221.0655 260.4775
240.7497 220.0558 261.4437
240.7309 219.0875 262.3743
240.7214 218.1604 263.2823
240.7083 217.2595 264.1571
240.7012 216.3942 265.0081
Code 
forecasts_data_xts_y10 <- TMS.arimapart[281:300, "y10"]
forecasts_xts_2 <- xts(forecasts_data,
                       order.by = index(forecasts_data_xts_y10))
forecasts_xts_2
             f_mean  f_lower  f_upper
2022-03-08 243.2799 241.9221 244.6377
2022-03-09 242.7846 240.1121 245.4571
2022-03-10 242.3408 238.0799 246.6016
2022-03-11 241.9608 236.2004 247.7211
2022-03-12 241.7312 234.4580 249.0043
2022-03-13 241.4810 232.7103 250.2517
2022-03-14 241.3199 231.1416 251.4983
2022-03-15 241.1933 229.6292 252.7574
2022-03-16 241.0672 228.1877 253.9467
2022-03-17 240.9990 226.8621 255.1359
2022-03-18 240.9235 225.5713 256.2756
2022-03-19 240.8691 224.3640 257.3742
2022-03-20 240.8340 223.2169 258.4511
2022-03-21 240.7928 222.1090 259.4765
2022-03-22 240.7715 221.0655 260.4775
2022-03-23 240.7497 220.0558 261.4437
2022-03-24 240.7309 219.0875 262.3743
2022-03-25 240.7214 218.1604 263.2823
2022-03-26 240.7083 217.2595 264.1571
2022-03-27 240.7012 216.3942 265.0081
Code 
forecasting_y10 <- merge(financial_data$y10, forecasts_xts_2)
tail(forecasting_y10, n = 20)
              y10   f_mean  f_lower  f_upper
2022-03-08 253.33 243.2799 241.9221 244.6377
2022-03-09 251.60 242.7846 240.1121 245.4571
2022-03-10 250.84 242.3408 238.0799 246.6016
2022-03-11 250.14 241.9608 236.2004 247.7211
2022-03-12 249.92 241.7312 234.4580 249.0043
2022-03-13 250.31 241.4810 232.7103 250.2517
2022-03-14 249.75 241.3199 231.1416 251.4983
2022-03-15 249.04 241.1933 229.6292 252.7574
2022-03-16 247.09 241.0672 228.1877 253.9467
2022-03-17 246.74 240.9990 226.8621 255.1359
2022-03-18 246.23 240.9235 225.5713 256.2756
2022-03-19 246.42 240.8691 224.3640 257.3742
2022-03-20 246.38 240.8340 223.2169 258.4511
2022-03-21 246.68 240.7928 222.1090 259.4765
2022-03-22 246.30 240.7715 221.0655 260.4775
2022-03-23 245.85 240.7497 220.0558 261.4437
2022-03-24 245.64 240.7309 219.0875 262.3743
2022-03-25 245.77 240.7214 218.1604 263.2823
2022-03-26 244.87 240.7083 217.2595 264.1571
2022-03-27 244.08 240.7012 216.3942 265.0081

Here we applied the same procedure as we did for y3.

Plotting forecast of y10

Code 
par(mfrow = c(1, 1))

plot(forecasting_y10[, c("y10", "f_mean", "f_lower", "f_upper")], 
     major.ticks = "days", 
     grid.ticks.on = "days",
     grid.ticks.lty = 3,
     main = "20 day forecast of y3",
     col = c("black", "blue", "red", "red"))

Code 
plot(forecasting_y10["2022-02/", # limit the rows
                    c("y10", "f_mean", "f_lower", "f_upper")], 
     major.ticks = "days", 
     grid.ticks.on = "days",
     grid.ticks.lty = 3,
     main = "20 day forecast of y10",
     col = c("black", "blue", "red", "red"))

Here we can visually interprete that previous ARIMA model was better y3, even though in the end of sample forecasting line approachs original line. However, we will make decision after comparing error metrics.

Error metrics for y3

Code 
TSy3 <- forecasting_y3
TSy3$mae <- abs(TSy3$y3 - TSy3$f_mean )
TSy3$mse <- (TSy3$y3 - TSy3$f_mean ) ^ 2
TSy3$mape <- abs((TSy3$y3 - TSy3$f_mean)/TSy3$y3)
TSy3$amape <-  abs((TSy3$y3 - TSy3$f_mean)/(TSy3$y3 + TSy3$f_mean))
tail(TSy3, 20)
               y3   f_mean  f_lower  f_upper       mae        mse        mape
2022-03-08 302.89 285.0073 283.0169 286.9976 17.882703 319.791074 0.059040256
2022-03-09 300.11 284.9394 280.3470 289.5318 15.170620 230.147709 0.050550198
2022-03-10 298.22 284.7407 277.1153 292.3662 13.479258 181.690390 0.045199040
2022-03-11 297.11 284.5761 273.9515 295.2007 12.533876 157.098036 0.042185977
2022-03-12 296.66 284.4397 270.8429 298.0365 12.220300 149.335727 0.041192947
2022-03-13 296.83 284.3266 267.8055 300.8478 12.503359 156.333989 0.042122963
2022-03-14 296.04 284.2329 264.8529 303.6130 11.807055 139.406552 0.039883310
2022-03-15 293.82 284.1553 261.9937 306.3169  9.664704  93.406507 0.032893282
2022-03-16 291.49 284.0909 259.2320 308.9498  7.399055  54.746008 0.025383562
2022-03-17 290.05 284.0376 256.5688 311.5065  6.012384  36.148759 0.020728784
2022-03-18 288.80 283.9934 254.0024 313.9844  4.806580  23.103207 0.016643281
2022-03-19 288.38 283.9568 251.5301 316.3835  4.423206  19.564752 0.015338117
2022-03-20 288.44 283.9264 249.1480 318.7049  4.513560  20.372221 0.015648175
2022-03-21 288.86 283.9013 246.8518 320.9508  4.958715  24.588851 0.017166498
2022-03-22 289.20 283.8804 244.6368 323.1240  5.319561  28.297734 0.018394058
2022-03-23 288.83 283.8632 242.4985 325.2278  4.966838  24.669479 0.017196406
2022-03-24 287.98 283.8488 240.4323 327.2654  4.131156  17.066446 0.014345286
2022-03-25 287.54 283.8370 238.4336 329.2404  3.703021  13.712364 0.012878281
2022-03-26 286.01 283.8271 236.4983 331.1560  2.182854   4.764853 0.007632091
2022-03-27 285.43 283.8190 234.6223 333.0157  1.611003   2.595332 0.005644128
                 amape
2022-03-08 0.030418074
2022-03-09 0.025930495
2022-03-10 0.023122068
2022-03-11 0.021547489
2022-03-12 0.021029610
2022-03-13 0.021514611
2022-03-14 0.020347416
2022-03-15 0.016721656
2022-03-16 0.012854933
2022-03-17 0.010472938
2022-03-18 0.008391471
2022-03-19 0.007728327
2022-03-20 0.007885787
2022-03-21 0.008657559
2022-03-22 0.009282399
2022-03-23 0.008672773
2022-03-24 0.007224462
2022-03-25 0.006480872
2022-03-26 0.003830663
2022-03-27 0.002830051
Code 
mean_errors <- colMeans(TSy3[, c("mae", "mse", "mape", "amape")], na.rm = TRUE)
print(mean_errors)
        mae         mse        mape       amape 
 7.96449033 84.84199955  0.02700333  0.01374718

Error metrics for y10

Code 
TSy10 <- forecasting_y10
TSy10$mae <- abs(TSy10$y10 - TSy10$f_mean )
TSy10$mse <- (TSy10$y10 - TSy10$f_mean ) ^ 2
TSy10$mape <- abs((TSy10$y10 - TSy10$f_mean)/TSy10$y10)
TSy10$amape <-  abs((TSy10$y10 - TSy10$f_mean)/(TSy10$y10 + TSy10$f_mean))
tail(TSy10, 20)
              y10   f_mean  f_lower  f_upper       mae       mse       mape
2022-03-08 253.33 243.2799 241.9221 244.6377 10.050080 101.00411 0.03967189
2022-03-09 251.60 242.7846 240.1121 245.4571  8.815409  77.71144 0.03503740
2022-03-10 250.84 242.3408 238.0799 246.6016  8.499246  72.23718 0.03388313
2022-03-11 250.14 241.9608 236.2004 247.7211  8.179237  66.89992 0.03269864
2022-03-12 249.92 241.7312 234.4580 249.0043  8.188822  67.05681 0.03276577
2022-03-13 250.31 241.4810 232.7103 250.2517  8.828979  77.95086 0.03527218
2022-03-14 249.75 241.3199 231.1416 251.4983  8.430059  71.06589 0.03375399
2022-03-15 249.04 241.1933 229.6292 252.7574  7.846698  61.57067 0.03150778
2022-03-16 247.09 241.0672 228.1877 253.9467  6.022838  36.27458 0.02437508
2022-03-17 246.74 240.9990 226.8621 255.1359  5.740982  32.95888 0.02326734
2022-03-18 246.23 240.9235 225.5713 256.2756  5.306546  28.15943 0.02155117
2022-03-19 246.42 240.8691 224.3640 257.3742  5.550929  30.81282 0.02252629
2022-03-20 246.38 240.8340 223.2169 258.4511  5.545964  30.75772 0.02250980
2022-03-21 246.68 240.7928 222.1090 259.4765  5.887248  34.65969 0.02386593
2022-03-22 246.30 240.7715 221.0655 260.4775  5.528499  30.56431 0.02244620
2022-03-23 245.85 240.7497 220.0558 261.4437  5.100262  26.01268 0.02074542
2022-03-24 245.64 240.7309 219.0875 262.3743  4.909076  24.09902 0.01998484
2022-03-25 245.77 240.7214 218.1604 263.2823  5.048623  25.48860 0.02054206
2022-03-26 244.87 240.7083 217.2595 264.1571  4.161712  17.31985 0.01699560
2022-03-27 244.08 240.7012 216.3942 265.0081  3.378820  11.41642 0.01384308
                 amape
2022-03-08 0.020237373
2022-03-09 0.017831076
2022-03-10 0.017233530
2022-03-11 0.016621062
2022-03-12 0.016655756
2022-03-13 0.017952704
2022-03-14 0.017166717
2022-03-15 0.016006049
2022-03-16 0.012337908
2022-03-17 0.011770603
2022-03-18 0.010892965
2022-03-19 0.011391451
2022-03-20 0.011383014
2022-03-21 0.012077081
2022-03-22 0.011350489
2022-03-23 0.010481433
2022-03-24 0.010093275
2022-03-25 0.010377621
2022-03-26 0.008570630
2022-03-27 0.006969783
Code 
mean_errors <- colMeans(TSy10[, c("mae", "mse", "mape", "amape")], na.rm = TRUE)
print(mean_errors)
        mae         mse        mape       amape 
 6.55100154 46.20104376  0.02636218  0.01337003

Overall, some error metrics(mape, amape) are almost the same. However, since lower values of other metrics (mae and mse) in forecasting y10, ARIMA model suits better. However, we cannot conclude that all error metrics are smaller than metrics of other model. I think main reason is that in y3 prices are higher y10 in testdata and and when we divide them as percentage mape and amape can show differences. If have to choose one of VECM and ARIMA models using error metrics, I would choose that VECM model suits better than ARIMA model