CoCalc Public Filesr-marima.ipynb
Author: Harald Schilly
Views : 174
Compute Environment: Ubuntu 18.04 (Deprecated)

# Marima R library in CoCalc

Multivariate ARIMA and ARIMA-X estimation using Spliid's algorithm (marima()) and simulation (marima.sim()).

https://cran.r-project.org/package=marima

In [1]:
R.version

_ platform x86_64-pc-linux-gnu arch x86_64 os linux-gnu system x86_64, linux-gnu status major 3 minor 4.4 year 2018 month 03 day 15 svn rev 74408 language R version.string R version 3.4.4 (2018-03-15) nickname Someone to Lean On
In [2]:
library(marima)
data(austr)
series<-t(austr)[,1:90]
# Define marima model
Model5 <- define.model(kvar=7,ar=1,ma=1,rem.var=1,reg.var=6:7)
# Estimate marima model
Marima5 <- marima(series,Model5$ar.pattern,Model5$ma.pattern,penalty=1)
# Calculate residuals by filtering
Resid <- arma.filter(series, Marima5$ar.estimates, Marima5$ma.estimates)
# Compare residuals
plot(Marima5$residuals[2, 5:90], Resid$residuals[2, 5:90], xlab='marima residuals', ylab='arma.filter residuals')

All cases in data, 1 to 90 accepted for completeness. 90 7 = MARIMA - dimension of data arma.filter is being called indicators for means= 1 1 1 1 1 1 1 dim(yseries) 7 94
In [3]:
library(marima)
data(austr)
series<-austr
Model5 <- define.model(kvar=7, ar=1, ma=1, rem.var=1, reg.var=6:7)
Marima5 <- marima(ts(series[1:90, ]), Model5$ar.pattern, Model5$ma.pattern,
penalty=1)
nstart <- 90
nstep <- 10
cat("Calling arma.forecast.\n")
cat("In the example the input series is dim(length,kvar).\n")
cat("and of type ts() (timeseries) for illustration. \n")
Forecasts <- arma.forecast(series=ts(series), marima=Marima5,
nstart=nstart, nstep=nstep )
Year<-series[91:100,1]
One.step <- Forecasts$forecasts[, (nstart+1)] One.step Predict <- Forecasts$forecasts[ 2, 91:100]
Predict
stdv<-sqrt(Forecasts$pred.var[2, 2, ]) upper.lim=Predict+stdv*1.645 lower.lim=Predict-stdv*1.645 Out<-rbind(Year, Predict, upper.lim, lower.lim) print(Out) # plot results: plot(series[1:100, 1], Forecasts$forecasts[2, ], type='l', xlab='Year',
ylab='Rate of armed suicides', main='Prediction of suicides by firearms',
ylim=c(0.0, 4.1))
lines(series[1:90, 1], series[1:90, 2], type='p')
grid(lty=2, lwd=1, col='black')
Years<-2005:2014
lines(Years, Predict, type='l')
lines(Years, upper.lim, type='l')
lines(Years, lower.lim, type='l')
lines(c(2004.5, 2004.5), c(0.0, 2.0), lty = 2)

Input data is transformed from type = 'time series' (is.ts) to matrix).Sampling information is ignored. One line is one sample at a certain time point.All cases in data, 1 to 90 accepted for completeness. 90 7 = MARIMA - dimension of data Calling arma.forecast. In the example the input series is dim(length,kvar). and of type ts() (timeseries) for illustration. Input series is type ts (timeseries). It will be changed to matrix(...) using as.matrix(...). Data matrix is to be transposed. Done! 7 Variables in series , with required length = 100 . print AR-model: , , 1 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1 0 0 0 0 0 0 [2,] 0 1 0 0 0 0 0 [3,] 0 0 1 0 0 0 0 [4,] 0 0 0 1 0 0 0 [5,] 0 0 0 0 1 0 0 [6,] 0 0 0 0 0 1 0 [7,] 0 0 0 0 0 0 1 , , 2 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0 0.00000000 0.0000000 0.00000000 0.00000000 0.00000000 0.0000000 [2,] 0 -0.56193313 -0.8153282 0.03416575 0.00000000 0.57325886 0.0000000 [3,] 0 -0.06075435 -0.3170937 0.00000000 -0.06694902 0.09659622 0.0000000 [4,] 0 0.65334775 -1.6631316 -0.82681628 -0.47195995 -0.96665683 0.3252574 [5,] 0 0.00000000 0.0000000 0.00000000 -0.98009562 0.00000000 0.0000000 [6,] 0 0.00000000 0.0000000 0.00000000 0.00000000 0.00000000 0.0000000 [7,] 0 0.00000000 0.0000000 0.00000000 0.00000000 0.00000000 0.0000000 print MA-model: , , 1 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1 0 0 0 0 0 0 [2,] 0 1 0 0 0 0 0 [3,] 0 0 1 0 0 0 0 [4,] 0 0 0 1 0 0 0 [5,] 0 0 0 0 1 0 0 [6,] 0 0 0 0 0 1 0 [7,] 0 0 0 0 0 0 1 , , 2 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0 0.0000000 0 0.00000000 0.0000000 0 0 [2,] 0 0.0000000 0 0.13025444 -0.2351108 0 0 [3,] 0 0.0000000 0 0.03905407 0.0000000 0 0 [4,] 0 0.5857279 0 0.00000000 0.0000000 0 0 [5,] 0 0.0000000 0 0.00000000 -0.7163302 0 0 [6,] 0 0.0000000 0 0.00000000 0.0000000 0 0 [7,] 0 0.0000000 0 0.00000000 0.0000000 0 0 variable no. 6 not random and regressor. Variable no. 6 : future values seem OK. variable no. 7 not random and regressor. Variable no. 7 : future values seem OK. Constant = 1959.5 1.192804 0.0959164 2.100926 0.02197434 0.1 0.5 Calculation of forecasting variances.
1. 1959.5
2. 1.00077878331063
3. 0.197060905030581
4. 8.22351229278395
5. 1.28373981144414
6. 0.1
7. 0.5
1. 1.00077878331063
2. 1.06162305060619
3. 1.15060432312317
4. 1.255825600993
5. 1.37462646550024
6. 1.5057582811964
7. 1.64816732488692
8. 1.80084140982766
9. 1.96281212350129
10. 2.13316980065769
[,1] [,2] [,3] [,4] [,5] Year 2005.0000000 2006.0000000 2007.0000000 2008.000000 2009.0000000 Predict 1.0007788 1.0616231 1.1506043 1.255826 1.3746265 upper.lim 1.4853282 1.7435373 1.9073152 2.044681 2.1785422 lower.lim 0.5162294 0.3797088 0.3938934 0.466970 0.5707108 [,6] [,7] [,8] [,9] [,10] Year 2010.000000 2011.0000000 2012.0000000 2013.000000 2014.000000 Predict 1.505758 1.6481673 1.8008414 1.962812 2.133170 upper.lim 2.317438 2.4642980 2.6198021 2.783731 2.955526 lower.lim 0.694079 0.8320367 0.9818807 1.141893 1.310814
In [4]:
# Generate Y=series with 4 variables for illustration:
set.seed(4711)
Y<-matrix(round(100*rnorm(40)+10), nrow=4)
# Example 1: use of difference parameter: If
difference=c(2, 1, 2, 1, 3, 12)
difference
# the variable 2 is differenced
# twice, and variable 3 is differenced once with lag=12.
# Example 2:
poly <- define.dif(series=Y, difference=c(2, 1, 3, 1, 3, 1))
poly
# Generates a (4-variate) polynomial differencing array (with a leading
# unity matrix corresponding to lag=0, and (in the example) differencing
# of variable 2 for lag 1 and variable 3 for lag 1 but twice. Afterwards
# the series Y is differenced accordingly. Results in poly$series and # poly$dif.poly .
# Example 3: Generation and application of multivariate differencing
# polynomial. Re-use the 4-variate time series and use the
# differencing polynomial (ar-form):
# var=1, dif=1, var=2, dif=6, and var=3 and 4, no differencing.
dif.y <-define.dif(Y, c(1, 1, 2, 6, 3, 0, 4, 0))
# Now dif.y contains the differenced series and the differencing
# polynomial. Print the generated polynomial in short form:
short.form(dif.y$dif.poly) # Specifying no differencing (3, 0 and 4, 0) may be omitted: dif.y <-define.dif(Y, c(1, 1, 2, 6)) dif.y # Example 4: y<-matrix(round(rnorm(1200)*100+50), nrow=6) library(marima) difference<-c(3, 2, 4, 0, 5, 0, 6, 7) matrix(difference, nrow=2) Y<-define.dif(y, difference=difference) round(rowMeans(Y$dif.series), 2)
round(Y$averages, 2)  1. 2 2. 1 3. 2 4. 1 5. 3 6. 12 arma.filter is being called indicators for means= 1 1 1 1 dim(yseries) 4 16$y.dif
 6 14 -37 81 -127 -57 124 -19 198 -5 -4 -30 12 68 -147 -197 -142 298 -191 224 -370 362 -227 -38 -147 168 113 -16 -6 187 -66 177
$y.lost  192 -51 145.6 -288 137.6 -175.6 -31 -86$dif.poly
1. 1
2. 0
3. 0
4. 0
5. 0
6. 1
7. 0
8. 0
9. 0
10. 0
11. 1
12. 0
13. 0
14. 0
15. 0
16. 1
17. 0
18. 0
19. 0
20. 0
21. 0
22. -1
23. 0
24. 0
25. 0
26. 0
27. -2
28. 0
29. 0
30. 0
31. 0
32. 0
33. 0
34. 0
35. 0
36. 0
37. 0
38. 0
39. 0
40. 0
41. 0
42. 0
43. 1
44. 0
45. 0
46. 0
47. 0
48. 0
$averages 1. 12.6 2. 1.4 3. -7.6 4. 29.3$dif.series
 192 -51 6 14 -37 81 -127 -57 124 -19 145.6 -288 198 -5 -4 -30 12 68 -147 -197 137.6 -175.6 -142 298 -191 224 -370 362 -227 -38 -31 -86 -147 168 113 -16 -6 187 -66 177
arma.filter is being called indicators for means= 1 1 1 1 dim(yseries) 4 28
1. 1
2. 0
3. 0
4. 0
5. 0
6. 1
7. 0
8. 0
9. 0
10. 0
11. 1
12. 0
13. 0
14. 0
15. 0
16. 1
17. 1
18. 0
19. 0
20. 0
21. 0
22. 2
23. 0
24. 0
25. 0
26. 0
27. 0
28. 0
29. 0
30. 0
31. 0
32. 0
33. -1
34. 0
35. 0
36. 0
37. 0
38. 1
39. 0
40. 0
41. 0
42. 0
43. -1
44. 0
45. 0
46. 0
47. 0
48. -1
49. -1
50. 0
51. 0
52. 0
53. 0
54. 0
55. 0
56. 0
57. 0
58. 0
59. 0
60. 0
61. 0
62. 0
63. 0
64. 0
65. 0
66. 0
67. 0
68. 0
69. 0
70. -1
71. 0
72. 0
73. 0
74. 0
75. 0
76. 0
77. 0
78. 0
79. 0
80. 0
81. 0
82. 0
83. 0
84. 0
85. 0
86. -2
87. 0
88. 0
89. 0
90. 0
91. 0
92. 0
93. 0
94. 0
95. 0
96. 0
97. 0
98. 0
99. 0
100. 0
101. 0
102. -1
103. 0
104. 0
105. 0
106. 0
107. 0
108. 0
109. 0
110. 0
111. 0
112. 0
arma.filter is being called indicators for means= 1 1 1 1 dim(yseries) 4 24
$y.dif  -208 70 181 -143 -117 239 -106 -298 -111 32 -52 -174 -6 187 -66 177$y.lost
 179.4 -243 57 8 -51 118 145.6 -142.4 55.6 50.6 46.6 16.6 130 92 -88 30 -43 108 -31 -86 -147 168 113 -16
$dif.poly 1. 1 2. 0 3. 0 4. 0 5. 0 6. 1 7. 0 8. 0 9. 0 10. 0 11. 1 12. 0 13. 0 14. 0 15. 0 16. 1 17. -1 18. 0 19. 0 20. 0 21. 0 22. 0 23. 0 24. 0 25. 0 26. 0 27. 0 28. 0 29. 0 30. 0 31. 0 32. 0 33. 0 34. 0 35. 0 36. 0 37. 0 38. 0 39. 0 40. 0 41. 0 42. 0 43. 0 44. 0 45. 0 46. 0 47. 0 48. 0 49. 0 50. 0 51. 0 52. 0 53. 0 54. 0 55. 0 56. 0 57. 0 58. 0 59. 0 60. 0 61. 0 62. 0 63. 0 64. 0 65. 0 66. 0 67. 0 68. 0 69. 0 70. 0 71. 0 72. 0 73. 0 74. 0 75. 0 76. 0 77. 0 78. 0 79. 0 80. 0 81. 0 82. 0 83. 0 84. 0 85. 0 86. 0 87. 0 88. 0 89. 0 90. 0 91. 0 92. 0 93. 0 94. 0 95. 0 96. 0 97. 0 98. 0 99. 0 100. 0 101. 0 102. -1 103. 0 104. 0 105. 0 106. 0 107. 0 108. 0 109. 0 110. 0 111. 0 112. 0$averages
1. 12.6
2. 1.4
3. -7.6
4. 29.3
\$dif.series
 179.4 -243 57 8 -51 118 -208 70 181 -143 145.6 -142.4 55.6 50.6 46.6 16.6 -117 239 -106 -298 130 92 -88 30 -43 108 -111 32 -52 -174 -31 -86 -147 168 113 -16 -6 187 -66 177
 3 4 5 6 2 0 0 7
arma.filter is being called indicators for means= 1 1 1 1 1 1 dim(yseries) 6 220
1. 44.7
2. 45.72
3. 4.12
4. 0.72
5. 0.25
6. 2.13
1. 45.14
2. 45.38
3. 51.83
4. 52.97
5. 48.83
6. 49.66
In [ ]: