Libraries
library(vars)
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: urca
## Loading required package: lmtest
library(nlme)
library(zoo)
library(readxl)
library(data.table)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(forecast)
##
## Attaching package: 'forecast'
## The following object is masked from 'package:nlme':
##
## getResponse
Data Screening
setwd("C:/Users/Gautam/OneDrive/Ruchil/565/R_script")
covid19_india <- read_excel("covid19_india.xlsx")
covid19_us <- read_excel("covid19_us.xlsx")
covid19_us$date <- as.Date(covid19_us$date)
covid19_india$date <- as.Date(covid19_india$date)
head(covid19_india)
head(covid19_us)
diff(range(covid19_us$date))
## Time difference of 172 days
diff(range(covid19_india$date))
## Time difference of 172 days
range(covid19_us$date)
## [1] "2020-01-22" "2020-07-12"
range(covid19_india$date)
## [1] "2020-01-22" "2020-07-12"
###plot for positive cases
par(lwd = 1, cex = 0.8, las = 1)
plot(cases ~ date,covid19_us,type = 'o', pch = 19,ylab='Positive cases',xlab='Months',
main = 'Covid 19 cases in US till 19th July')
par(lwd = 1, cex = 0.8, las = 1)
plot(cases ~ date,covid19_india,type = 'o', pch = 19,
ylab='Positive cases',xlab='Months',
main= 'Covid-19 cases in India till 19th July')
###plot for deaths
par(lwd = 1, cex = 0.8, las = 1)
plot(death ~ date,covid19_us,type = 'o', pch = 19,ylab='Deaths',xlab='Months',
main = 'Deaths due to Covid-19 in US till 19th July')
par(lwd = 1, cex = 0.8, las = 1)
plot(death ~ date,covid19_india,type = 'o', pch = 19,ylab='Deaths',xlab='Months',
main = 'Deaths due to Covid-19 in India till 19th July')
Time series
usz <- zoo(
covid19_us[,-1],
covid19_us$date
)
inz <- zoo(
covid19_india[,-1],
covid19_india$date
)
### setting the start date as
start <- as.Date('2020-03-16')
newus <- usz[time(usz) >= start]
newin <- inz[time(inz) >= start]
#Time series
usts <- ts(newus,start = 11,freq = 7)
ints <- ts(newin,start = 11,freq = 7)
plot(usts, main = 'Time Series plot for USA covid-19 cases and deaths', xlab = 'week')
plot(ints, main = 'Time Series plot for INDIA covid-19 cases and deaths', xlab = 'week')
SEASONAL TIME DECOMPOSITION
decom.uscases <- decompose(usts[,1])
decom.usdeath <- decompose(usts[,2])
decom.incases <- decompose(ints[,1])
decom.indeath <- decompose(ints[,2])
plot(decom.uscases)
plot(decom.usdeath)
plot(decom.incases)
plot(decom.indeath)
uscasesrand <- na.omit(decom.uscases$random)
usdeathrand <- na.omit(decom.usdeath$random)
incasesrand <- na.omit(decom.incases$random)
indeathrand <- na.omit(decom.indeath$random)
ccf(uscasesrand,incasesrand)
ccf(uscasesrand,usdeathrand)
ccf(incasesrand,indeathrand)
plot(aggregate(usts))
plot(aggregate(ints))
#for both india and US the direction of the commulative plot is increasing so the deaths and cases will increase
Type of Series - Stationary or non Stationary ?
casesus <- usts[,1]#non stationary
deathus <- usts[,2]#non stationary
casesin <- ints[,1]#non stationary
deathin <- ints[,2]#stationary
adf.test(casesus)$p.value
## Warning in adf.test(casesus): p-value greater than printed p-value
## [1] 0.99
adf.test(deathus)$p.value
## [1] 0.3339229
adf.test(casesin)$p.value
## Warning in adf.test(casesin): p-value greater than printed p-value
## [1] 0.99
adf.test(deathin)$p.value
## [1] 0.02538913
acf(casesus)
acf(diff(casesus))
pacf(casesus)
acf(deathus)
acf(diff(deathus))
pacf(deathus)
acf(casesin)
acf(diff(casesin))
pacf(casesin)
acf(deathin)
acf(diff(deathin))
pacf(deathin)
Making all series stationary
casesusd <- log(casesus)
deathusd <- diff((deathus))
casesind <- log((casesin))
plot(casesusd)
plot(deathusd)
plot(casesind)
acf(casesusd)
pacf(casesusd)
acf(deathusd)
pacf(deathusd)
acf(casesind)
pacf(casesind)
adf.test(casesusd)$p.value
## [1] 0.05067324
adf.test(deathusd)$p.value
## Warning in adf.test(deathusd): p-value smaller than printed p-value
## [1] 0.01
adf.test(casesind)$p.value
## [1] 0.01595341
US cases Analysis
casesusd <- log(casesus)
Time <- time(casesusd)
Seas<- factor(cycle(casesusd))
casesusd.lm <-lm(casesusd ~ 0+Time + I(Time^2)+I(Time^3)+I(Time^4)+Seas)
summary(casesusd.lm)
##
## Call:
## lm(formula = casesusd ~ 0 + Time + I(Time^2) + I(Time^3) + I(Time^4) +
## Seas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48240 -0.08575 0.00508 0.08511 0.34892
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Time 1.655e+01 8.338e-01 19.84 <2e-16 ***
## I(Time^2) -1.215e+00 6.761e-02 -17.98 <2e-16 ***
## I(Time^3) 3.842e-02 2.373e-03 16.19 <2e-16 ***
## I(Time^4) -4.411e-04 3.049e-05 -14.47 <2e-16 ***
## Seas1 -7.177e+01 3.749e+00 -19.14 <2e-16 ***
## Seas2 -7.169e+01 3.749e+00 -19.12 <2e-16 ***
## Seas3 -7.164e+01 3.750e+00 -19.11 <2e-16 ***
## Seas4 -7.155e+01 3.750e+00 -19.08 <2e-16 ***
## Seas5 -7.150e+01 3.750e+00 -19.07 <2e-16 ***
## Seas6 -7.155e+01 3.750e+00 -19.08 <2e-16 ***
## Seas7 -7.168e+01 3.750e+00 -19.12 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1321 on 108 degrees of freedom
## Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
## F-statistic: 6.343e+04 on 11 and 108 DF, p-value: < 2.2e-16
acf(resid(casesusd.lm))[1]
##
## Autocorrelations of series 'resid(casesusd.lm)', by lag
##
## 1
## 0.569
pacf(resid(casesusd.lm))
new.Time <- seq(28, len = 7, by = 1/7)
new.Seas <- c(1,2,3,4,5,6,7)
new.data <- data.frame(Time=new.Time,Seas=factor(new.Seas))
predus.lm <- exp(predict(casesusd.lm,new.data)[1:7])
predus.lm
## 1 2 3 4 5 6 7
## 54375.46 59912.09 63994.42 70792.08 74440.91 71028.99 62420.36
uscases.gls <- gls(casesusd ~Time+I(Time^2) + I(Time^3)+I(Time^4) +Seas, cor = corAR1(0.569))
summary(uscases.gls)
## Generalized least squares fit by REML
## Model: casesusd ~ Time + I(Time^2) + I(Time^3) + I(Time^4) + Seas
## Data: NULL
## AIC BIC logLik
## -96.49174 -61.62403 61.24587
##
## Correlation Structure: AR(1)
## Formula: ~1
## Parameter estimate(s):
## Phi
## 0.7659299
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -83.59770 9.848790 -8.488118 0.0000
## Time 19.03320 2.206047 8.627738 0.0000
## I(Time^2) -1.40580 0.179669 -7.824410 0.0000
## I(Time^3) 0.04472 0.006319 7.077404 0.0000
## I(Time^4) -0.00052 0.000081 -6.374579 0.0000
## Seas2 0.07744 0.026147 2.961668 0.0038
## Seas3 0.12681 0.033408 3.795774 0.0002
## Seas4 0.21426 0.036442 5.879497 0.0000
## Seas5 0.25414 0.036531 6.956771 0.0000
## Seas6 0.19999 0.033700 5.934418 0.0000
## Seas7 0.06673 0.026778 2.491920 0.0142
##
## Correlation:
## (Intr) Time I(T^2) I(T^3) I(T^4) Seas2 Seas3 Seas4 Seas5 Seas6
## Time -0.998
## I(Time^2) 0.991 -0.998
## I(Time^3) -0.981 0.992 -0.998
## I(Time^4) 0.969 -0.983 0.993 -0.998
## Seas2 0.048 -0.048 0.048 -0.047 0.047
## Seas3 0.060 -0.061 0.060 -0.059 0.058 0.639
## Seas4 0.064 -0.064 0.062 -0.061 0.060 0.470 0.723
## Seas5 0.060 -0.059 0.057 -0.055 0.053 0.361 0.547 0.743
## Seas6 0.046 -0.044 0.041 -0.039 0.036 0.275 0.410 0.549 0.726
## Seas7 0.014 -0.011 0.007 -0.004 0.000 0.189 0.279 0.367 0.478 0.648
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.88903944 -0.51777854 0.06121349 0.59275193 2.63195172
##
## Residual standard error: 0.1620561
## Degrees of freedom: 119 total; 108 residual
acf(resid(uscases.gls))[1]
##
## Autocorrelations of series 'resid(uscases.gls)', by lag
##
## 1
## 0.369
pacf(resid(uscases.gls))
predglsusc <- exp(predict(uscases.gls,new.data))
predglsusc
## [1] 54303.01 59089.25 62291.29 67952.59 70393.84 66089.95 57065.35
## attr(,"label")
## [1] "Predicted values"
AIC(uscases.gls)
## [1] -96.49174
AIC(casesusd.lm)
## [1] -131.5131
get.best.arima <- function(x.ts, maxord = c(1,1,1,1,1,1))
{
best.aic <- 1e8
n <- length(x.ts)
for (p in 0:maxord[1]) for(d in 0:maxord[2]) for(q in 0:maxord[3])
for (P in 0:maxord[4]) for(D in 0:maxord[5]) for(Q in 0:maxord[6])
{
fit <- arima(x.ts, order = c(p,d,q),
seas = list(order = c(P,D,Q),
frequency(x.ts)), method = "CSS")
fit.aic <- -2 * fit$loglik + (log(n) + 1) * length(fit$coef)
if (fit.aic < best.aic)
{
best.aic <- fit.aic
best.fit <- fit
best.model <- c(p,d,q,P,D,Q)
}
}
list(best.aic, best.fit, best.model)
}
us.cases.arima <- get.best.arima(resid(uscases.gls), maxord = c(2,2,2,2,2,2))
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
us.cases.arima
## [[1]]
## [1] -262.5335
##
## [[2]]
##
## Call:
## arima(x = x.ts, order = c(p, d, q), seasonal = list(order = c(P, D, Q), frequency(x.ts)),
## method = "CSS")
##
## Coefficients:
## ar1 ar2 sar1
## -0.4539 -0.2162 0.1250
## s.e. 0.0910 0.0827 0.0841
##
## sigma^2 estimated as 0.005464: part log likelihood = 139.94
##
## [[3]]
## [1] 2 1 0 1 0 0
us.cases.pred <-predict(us.cases.arima[[2]],n.ahead=7)
acf(us.cases.arima[[2]]$residuals)
pacf(us.cases.arima[[2]]$residuals)
us.cases.pred <- us.cases.pred$pred
us.cases.pred <- exp(us.cases.pred)
preduscaseslm <- ts(c(predus.lm+us.cases.pred),start=28,frequency=7)
preduscaseslm
## Time Series:
## Start = c(28, 1)
## End = c(28, 7)
## Frequency = 7
## 1 2 3 4 5 6 7
## 54376.50 59913.12 63995.48 70793.11 74441.95 71030.02 62421.41
actualuscases <- ts(c(58465,62879,65382,70953,77233,65180,63907),start=28,frequency=7)
ts.plot(actualuscases,preduscaseslm,col=c("red","blue"),lty = 1:2,ylab= 'Positive cases',main='Prediction of US covid-19 cases from 13th to 19th July 2020 ')
ts.plot(cbind(window(casesus, start = 10), preduscaseslm),col=c("red","blue"), lty = 1:2,ylab= 'Positive covid-19 cases',main='US covid-19 cases from 16th Mar to 19th July 2020 ',xlab='weeks')
## Warning in window.default(x, ...): 'start' value not changed
mpe <- function(actual,predict){
mpe <- mean(abs((actual - predict)/actual))*100
return (mpe)
}
mpe(actualuscases,preduscaseslm)
## [1] 4.138498
US deaths Analysis
deathusd <- diff(deathus)
casesusdi <- diff(casesus)
ustsnew <- ts(cbind(deathusd,casesusdi),start = 11,freq = 7)
usnewcc <- na.omit(ustsnew)
ddeath.us <- usnewcc[,1]
casesd.us <- usnewcc[,2]
#time(ddeath.us) last time = 27.71429 net will be 27.85714
po.test(cbind(ddeath.us,casesd.us))
## Warning in po.test(cbind(ddeath.us, casesd.us)): p-value smaller than printed p-
## value
##
## Phillips-Ouliaris Cointegration Test
##
## data: cbind(ddeath.us, casesd.us)
## Phillips-Ouliaris demeaned = -145.6, Truncation lag parameter = 1,
## p-value = 0.01
cases.deathus.vars <- VAR(cbind(ddeath.us, casesd.us), p = 3, type = "trend")
summary(cases.deathus.vars)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: ddeath.us, casesd.us
## Deterministic variables: trend
## Sample size: 115
## Log Likelihood: -1906.962
## Roots of the characteristic polynomial:
## 0.8378 0.8378 0.6692 0.6692 0.669 0.669
## Call:
## VAR(y = cbind(ddeath.us, casesd.us), p = 3, type = "trend")
##
##
## Estimation results for equation ddeath.us:
## ==========================================
## ddeath.us = ddeath.us.l1 + casesd.us.l1 + ddeath.us.l2 + casesd.us.l2 + ddeath.us.l3 + casesd.us.l3 + trend
##
## Estimate Std. Error t value Pr(>|t|)
## ddeath.us.l1 -0.34205 0.09200 -3.718 0.000320 ***
## casesd.us.l1 0.01371 0.01108 1.237 0.218592
## ddeath.us.l2 -0.09128 0.09847 -0.927 0.356026
## casesd.us.l2 -0.03054 0.01079 -2.831 0.005528 **
## ddeath.us.l3 -0.15649 0.09546 -1.639 0.104038
## casesd.us.l3 -0.04172 0.01065 -3.916 0.000158 ***
## trend 0.28882 0.49240 0.587 0.558724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 348.4 on 108 degrees of freedom
## Multiple R-Squared: 0.3019, Adjusted R-squared: 0.2567
## F-statistic: 6.674 on 7 and 108 DF, p-value: 1.445e-06
##
##
## Estimation results for equation casesd.us:
## ==========================================
## casesd.us = ddeath.us.l1 + casesd.us.l1 + ddeath.us.l2 + casesd.us.l2 + ddeath.us.l3 + casesd.us.l3 + trend
##
## Estimate Std. Error t value Pr(>|t|)
## ddeath.us.l1 3.18607 0.78426 4.063 9.23e-05 ***
## casesd.us.l1 -0.13895 0.09442 -1.472 0.144003
## ddeath.us.l2 3.18212 0.83939 3.791 0.000248 ***
## casesd.us.l2 -0.23502 0.09195 -2.556 0.011977 *
## ddeath.us.l3 2.42945 0.81370 2.986 0.003501 **
## casesd.us.l3 -0.25129 0.09081 -2.767 0.006652 **
## trend 12.45755 4.19744 2.968 0.003693 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 2970 on 108 degrees of freedom
## Multiple R-Squared: 0.3193, Adjusted R-squared: 0.2752
## F-statistic: 7.236 on 7 and 108 DF, p-value: 4.238e-07
##
##
##
## Covariance matrix of residuals:
## ddeath.us casesd.us
## ddeath.us 120921 294787
## casesd.us 294787 8819136
##
## Correlation matrix of residuals:
## ddeath.us casesd.us
## ddeath.us 1.0000 0.2855
## casesd.us 0.2855 1.0000
acf(cases.deathus.vars$varresult$ddeath.us$residuals)
acf(cases.deathus.vars$varresult$casesd.us$residuals)
cases.deathus.vars.pred <- predict(cases.deathus.vars, n.ahead = 9)
ddeath.us.pred <- ts(cases.deathus.vars.pred$fcst$ddeath.us[, 1], st = 27.85714, fr = 7)
casesd.us.pred <- ts(cases.deathus.vars.pred$fcst$casesd.us[, 1], st = 27.85714, fr = 7)
AIC(cases.deathus.vars)
## [1] 3841.925
uscasesprednew <- ts(c(
63007-578.47924 , 63007-578.47924+1514.30655,63007-578.47924+1514.30655+ 1944.75540 ,63007-578.47924+1514.30655+1944.75540+2127.31979,63007-578.47924+1514.30655+1944.75540+2127.31979+ 1467.71600,63007-578.47924+1514.30655+1944.75540+2127.31979+ 1467.71600 +213.01904 , 63007-578.47924+1514.30655+1944.75540+2127.31979+ 1467.71600+213.01904+ 52.94763,463.59088+63007-578.47924+1514.30655+1944.75540+2127.31979+ 1467.71600+213.01904+ 52.94763),start=27.85714,fr=7)
actualuscases12 <- ts(c(60978
,58465,62879,65382,70953,77233,65180,63907),start=27.85714,frequency=7)
ts.plot(actualuscases12,uscasesprednew,col=c("red","blue"),lty = 1:2,ylab= 'Positive cases',main='US covid-19 cases from 12 July to 19th July 2020 by vars() ')
usdeathprednew <- ts(c(374.8731,690.7686,793.9434,785.4235,671.5968, 604.9877,544.9384),start=28,frequency=7)
actualusdeath12 <- ts(c(327,736,1021,979,966,823,420),start=28,frequency=7)
ts.plot(actualusdeath12,usdeathprednew,col=c("red","blue"),lty = 1:2,ylab= 'Deaths',main='US deaths from 13 July to 19th July 2020 ')
ts.plot(cbind(deathus, usdeathprednew),col=c("red","blue"), lty = 1:2,ylab= 'Deaths',main='US covid-19 cases from 16th Mar to 19th July 2020 ',xlab='Week')
mpe(actualusdeath12,usdeathprednew)
## [1] 21.3587
India cases Analysis
casesind <- log((casesin))
Time <- time(casesind)
Seas<- factor(cycle(casesind))
casesind.lm <-lm(casesind ~ 0+Time + I(Time^2)+I(Time^3)+I(Time^4)+Seas)
summary(casesind.lm)
##
## Call:
## lm(formula = casesind ~ 0 + Time + I(Time^2) + I(Time^3) + I(Time^4) +
## Seas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71322 -0.08552 -0.01037 0.06747 0.76200
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Time 1.034e+01 1.410e+00 7.332 4.39e-11 ***
## I(Time^2) -6.874e-01 1.144e-01 -6.010 2.55e-08 ***
## I(Time^3) 2.074e-02 4.015e-03 5.165 1.11e-06 ***
## I(Time^4) -2.340e-04 5.158e-05 -4.536 1.49e-05 ***
## Seas1 -5.183e+01 6.342e+00 -8.172 6.25e-13 ***
## Seas2 -5.189e+01 6.343e+00 -8.182 5.96e-13 ***
## Seas3 -5.181e+01 6.343e+00 -8.168 6.38e-13 ***
## Seas4 -5.185e+01 6.343e+00 -8.173 6.22e-13 ***
## Seas5 -5.176e+01 6.343e+00 -8.160 6.68e-13 ***
## Seas6 -5.173e+01 6.343e+00 -8.155 6.82e-13 ***
## Seas7 -5.177e+01 6.343e+00 -8.161 6.63e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2235 on 108 degrees of freedom
## Multiple R-squared: 0.9993, Adjusted R-squared: 0.9992
## F-statistic: 1.428e+04 on 11 and 108 DF, p-value: < 2.2e-16
acf(resid(casesind.lm))[1]
##
## Autocorrelations of series 'resid(casesind.lm)', by lag
##
## 1
## 0.133
pacf(resid(casesind.lm))
new.Time <- seq(28, len = 7, by = 1/7)
new.Seas <- c(1,2,3,4,5,6,7)
new.data <- data.frame(Time=new.Time,Seas=factor(new.Seas))
predin.lm <- exp(predict(casesind.lm,new.data))
predin.lm
## 1 2 3 4 5 6 7
## 25915.12 24505.85 26777.21 25993.59 28428.07 29218.37 28173.88
AIC(casesind.lm)
## [1] -6.394846
in.cases.arima <- get.best.arima(resid(casesind.lm), maxord = c(2,2,2,2,2,2))
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
## Warning in arima(x.ts, order = c(p, d, q), seas = list(order = c(P, D, Q), :
## possible convergence problem: optim gave code = 1
in.cases.arima
## [[1]]
## [1] -51.90654
##
## [[2]]
##
## Call:
## arima(x = x.ts, order = c(p, d, q), seasonal = list(order = c(P, D, Q), frequency(x.ts)),
## method = "CSS")
##
## Coefficients:
## sar1 sar2 sma1
## 0.1978 0.1633 -1.0822
## s.e. 0.0008 0.0008 0.0013
##
## sigma^2 estimated as 0.03256: part log likelihood = 34.62
##
## [[3]]
## [1] 0 0 0 2 1 1
in.cases.pred <-predict(in.cases.arima[[2]],n.ahead=7)
## Warning in predict.Arima(in.cases.arima[[2]], n.ahead = 7): seasonal MA part of
## model is not invertible
acf(in.cases.arima[[2]]$residuals)
pacf(in.cases.arima[[2]]$residuals)
in.cases.pred <- in.cases.pred$pred
in.cases.pred <- exp(in.cases.pred)
predincaseslm <- ts(c(predin.lm+in.cases.pred),start=28,frequency=7)
predincaseslm
## Time Series:
## Start = c(28, 1)
## End = c(28, 7)
## Frequency = 7
## 1 2 3 4 5 6 7
## 25916.14 24506.87 26778.22 25994.60 28429.08 29219.38 28174.89
actualincases <- ts(c(28498,29429,32695,34956,34884,38902,40425),start=28,frequency=7)
ts.plot(actualincases,predincaseslm,col=c("red","blue"),lty = 1:2,ylab= 'Positive cases',main='INDIA covid-19 cases from 13th to 19th July 2020 ')
ts.plot(cbind(window(casesin, start = 11), predincaseslm),col=c("red","blue"), lty = 1:2,ylab= 'Positive covid-19 cases',main='INDIA covid-19 cases from 16th Mar to 19th July 2020 ',xlab='weeks')
mpe(actualincases,predincaseslm)
## [1] 20.45934
India Deaths Analysis
deathind <- diff(deathin)
casesindi <- diff(casesin)
intsnew <- ts(cbind(deathind,casesindi),start = 11,freq = 7)
innewcc <- na.omit(intsnew)
ddeath.in <- innewcc[,1]
casesd.in <- innewcc[,2]
#time(ddeath.in) last time = 27.71429 net will be 27.85714
po.test(cbind(ddeath.in,casesd.in))
## Warning in po.test(cbind(ddeath.in, casesd.in)): p-value smaller than printed p-
## value
##
## Phillips-Ouliaris Cointegration Test
##
## data: cbind(ddeath.in, casesd.in)
## Phillips-Ouliaris demeaned = -166.73, Truncation lag parameter = 1,
## p-value = 0.01
cases.deathin.vars <- VAR(cbind(ddeath.in, casesd.in), p = 3, type = "trend")
summary(cases.deathin.vars)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: ddeath.in, casesd.in
## Deterministic variables: trend
## Sample size: 115
## Log Likelihood: -1639.039
## Roots of the characteristic polynomial:
## 0.7257 0.7257 0.6709 0.6709 0.6159 0.5555
## Call:
## VAR(y = cbind(ddeath.in, casesd.in), p = 3, type = "trend")
##
##
## Estimation results for equation ddeath.in:
## ==========================================
## ddeath.in = ddeath.in.l1 + casesd.in.l1 + ddeath.in.l2 + casesd.in.l2 + ddeath.in.l3 + casesd.in.l3 + trend
##
## Estimate Std. Error t value Pr(>|t|)
## ddeath.in.l1 -0.778861 0.092888 -8.385 2.10e-13 ***
## casesd.in.l1 -0.070077 0.029437 -2.381 0.01904 *
## ddeath.in.l2 -0.511392 0.108378 -4.719 7.15e-06 ***
## casesd.in.l2 -0.030712 0.030307 -1.013 0.31314
## ddeath.in.l3 -0.261400 0.093808 -2.787 0.00629 **
## casesd.in.l3 0.000391 0.030790 0.013 0.98989
## trend 0.633450 0.325581 1.946 0.05430 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 175.8 on 108 degrees of freedom
## Multiple R-Squared: 0.4125, Adjusted R-squared: 0.3745
## F-statistic: 10.83 on 7 and 108 DF, p-value: 2.727e-10
##
##
## Estimation results for equation casesd.in:
## ==========================================
## casesd.in = ddeath.in.l1 + casesd.in.l1 + ddeath.in.l2 + casesd.in.l2 + ddeath.in.l3 + casesd.in.l3 + trend
##
## Estimate Std. Error t value Pr(>|t|)
## ddeath.in.l1 0.21446 0.29090 0.737 0.46259
## casesd.in.l1 0.14137 0.09219 1.534 0.12807
## ddeath.in.l2 -0.38830 0.33941 -1.144 0.25513
## casesd.in.l2 -0.18086 0.09491 -1.906 0.05937 .
## ddeath.in.l3 -0.21966 0.29378 -0.748 0.45626
## casesd.in.l3 -0.30992 0.09643 -3.214 0.00173 **
## trend 6.13651 1.01964 6.018 2.46e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 550.5 on 108 degrees of freedom
## Multiple R-Squared: 0.3686, Adjusted R-squared: 0.3277
## F-statistic: 9.006 on 7 and 108 DF, p-value: 1.027e-08
##
##
##
## Covariance matrix of residuals:
## ddeath.in casesd.in
## ddeath.in 30891 -7076
## casesd.in -7076 301249
##
## Correlation matrix of residuals:
## ddeath.in casesd.in
## ddeath.in 1.00000 -0.07335
## casesd.in -0.07335 1.00000
AIC(cases.deathin.vars)
## [1] 3306.078
acf(cases.deathin.vars$varresult$ddeath.in$residuals)
acf(cases.deathin.vars$varresult$casesd.in$residuals)
cases.deathin.vars.pred <- predict(cases.deathin.vars,8)
ddeath.in.pred <- ts(cases.deathin.vars.pred$fcst$ddeath.in[, 1], st = 27.85714, fr = 7)
casesd.in.pred <- ts(cases.deathin.vars.pred$fcst$casesd.in[, 1], st = 27.85714, fr = 7)
incasesprednew <- ts(c(29108+299.9249+236.0146,29108+299.9249+236.0146+595.2179,29108+299.9249+236.0146+679.0945+595.2179,29108+299.9249+236.0146+679.0945+595.2179+647.5364,29108+299.9249+236.0146+679.0945+595.2179+647.5364+547.1482,29108+299.9249+236.0146+679.0945+595.2179+647.5364+547.148+523.2495,29108+299.9249+236.0146+679.0945+595.2179+647.5364+547.148+523.2495+544.4187,29108+299.9249+236.0146+679.0945+595.2179+647.5364+547.148+523.2495+544.418+585.4392),start=27.85714,fr=7)
actualincases12 <- ts(c(28701,28498,29429,32695,34956,34884,38902,40425
),start=27.85714,frequency=7)
ts.plot(actualincases12,incasesprednew,col=c("red","blue"),lty = 1:2,ylab= 'Positive cases',main='INDIA covid-19 cases from 12 July to 19th July 2020 by vars() ')
indeathprednew <- ts(c( 551+7.6929392+ 52.8506576 +19.3150536,551+7.6929392+ 52.8506576 -15.6434041+19.3150536, 551+7.6929392+ 52.8506576 -15.6434041+19.3150536+0.6286506,551+7.6929392+ 52.8506576 -15.6434041+19.3150536+0.6286506+15.0080100,551+7.6929392+ 52.8506576 -15.6434041+19.3150536+0.6286506+15.0080100+13.2955839,551+7.6929392+ 52.8506576 -15.6434041+19.3150536+0.6286506+15.0080100+13.2955839 +8.401331,551+7.6929392+ 52.8506576 -15.6434041+19.3150536+0.6286506+15.0080100+13.2955839 +8.401331+9.1749420),start=28,frequency=7)
actualindeath12 <- ts(c(553,582,606,687,671,543,681),start=28,frequency=7)
ts.plot(actualindeath12,indeathprednew,col=c("red","blue"),lty = 1:2,ylab= 'Deaths',main='INIDA covid-19 deaths from 13th to 19th July 2020 ')
ts.plot(cbind(window(deathin, start = 11), indeathprednew),col=c("red","blue"), lty = 1:2,ylab= 'Deaths',main='INDIA covid-19 deaths from 16th Mar to 19th July 2020 ',xlab='weeks')
mpe(actualindeath12,indeathprednew)
## [1] 8.084419
Deaths in India Analysis (Another Method)
Time <- time(deathin)
Seas = factor(cycle(deathin))
length(casesin)
## [1] 119
length(Time)
## [1] 119
length(deathin)
## [1] 119
deathin.lmts <- lm(deathin ~ 0+Time+Seas)
summary(deathin.lmts)
##
## Call:
## lm(formula = deathin ~ 0 + Time + Seas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -177.78 -56.86 -15.46 28.38 1549.13
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Time 35.362 3.031 11.666 < 2e-16 ***
## Seas1 -505.709 69.717 -7.254 5.82e-11 ***
## Seas2 -399.878 70.075 -5.706 9.72e-08 ***
## Seas3 -504.342 70.434 -7.161 9.27e-11 ***
## Seas4 -512.394 70.794 -7.238 6.30e-11 ***
## Seas5 -502.387 71.154 -7.061 1.53e-10 ***
## Seas6 -504.556 71.516 -7.055 1.57e-10 ***
## Seas7 -516.961 71.878 -7.192 7.92e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 162 on 111 degrees of freedom
## Multiple R-squared: 0.7385, Adjusted R-squared: 0.7197
## F-statistic: 39.19 on 8 and 111 DF, p-value: < 2.2e-16
AIC(deathin.lmts)
## [1] 1558.259
acf(resid(deathin.lmts))[1]
##
## Autocorrelations of series 'resid(deathin.lmts)', by lag
##
## 1
## 0.088
pacf(resid(deathin.lmts))
new.timein <- seq(28, len = 7, by = 1/7)
alpha2 <- coef(deathin.lmts)[1]
beta1 <- coef(deathin.lmts)[2]
beta2 <- coef(deathin.lmts)[3]
beta3 <- coef(deathin.lmts)[4]
beta4 <- coef(deathin.lmts)[5]
beta5 <- coef(deathin.lmts)[6]
beta6 <- coef(deathin.lmts)[7]
beta7 <- coef(deathin.lmts)[8]
forc20191 <- alpha2*new.Time[1] + beta1
forc20192 <- alpha2*new.Time[2] + beta2
forc20193 <- alpha2*new.Time[3] + beta3
forc20194 <- alpha2*new.Time[4] + beta4
forc20195 <- alpha2*new.Time[5] + beta5
forc20196 <- alpha2*new.Time[6] + beta6
forc20197 <- alpha2*new.Time[7] + beta7
indeathforc <- ts(c(forc20191,forc20192,forc20193,forc20194,forc20195,forc20196,forc20197),start=28,fr=7)
actualindeath12 <- ts(c(553,582,606,687,671,543,681),start=28,frequency=7)
ts.plot(actualindeath12,indeathforc,col=c("red","blue"),lty = 1:2,ylab= 'Deaths',main='INIDA covid-19 deaths from 13th to 19th July 2020 ')
ts.plot(cbind(window(deathin, start = 11), indeathforc),col=c("red","blue"), lty = 1:2,ylab= 'Deaths',main='INDIA covid-19 deaths from 16th Mar to 19th July 2020 ',xlab='weeks')
mpe(actualindeath12,indeathforc)
## [1] 16.76982
casesin.ar <- ar(casesin, FUN = mean, method = "mle")
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
casesin.ar <- ar(casesin, method = "mle")
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
acf(na.omit(casesin.ar$resid))
pacf(na.omit(casesin.ar$resid))
predict(casesin.ar, n.ahead=7)
## $pred
## Time Series:
## Start = c(28, 1)
## End = c(28, 7)
## Frequency = 7
## [1] 28669.47 28552.77 29490.84 30972.09 32302.18 33167.72 33250.66
##
## $se
## Time Series:
## Start = c(28, 1)
## End = c(28, 7)
## Frequency = 7
## [1] 480.0993 683.6045 782.1472 839.3641 917.1031 989.5954 1118.9370
actualincases <- ts(c(28701,28498,29429,32695,34956,34884,38902),start=28,frequency=7)
Prediction of Deaths in India
deathin.ar <- ar(deathin, FUN = mean, method = "mle")
deathin.ar <- ar(deathin, method = "mle")
deathin.ar
##
## Call:
## ar(x = deathin, method = "mle")
##
## Coefficients:
## 1 2 3 4 5
## 0.1860 0.1511 0.1615 0.1609 0.2290
##
## Order selected 5 sigma^2 estimated as 28026
acf(na.omit(deathin.ar$resid))
pacf(na.omit(deathin.ar$resid))
prednew<- predict(deathin.ar, n.ahead=7)
actualindeath12 <- ts(c(553,582,606,687,671,543,681),start=28,frequency=7)
mpe(actualindeath12,prednew$pred)
## [1] 25.59061
Prediction of Deaths in India
model.arma <- auto.arima(deathin, ic="aic")
model.arma
## Series: deathin
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## -0.9113 4.6783
## s.e. 0.0429 1.4563
##
## sigma^2 estimated as 26408: log likelihood=-768.02
## AIC=1542.03 AICc=1542.24 BIC=1550.35
acf(resid(model.arma))
pacf(resid(model.arma))
new.time <- seq(length(deathin), length = 7)
new.data <- data.frame(Time = new.time, Imth = 1:7)
arma<-arima(deathin, order=c(0,1,1), include.mean = FALSE)
arma
##
## Call:
## arima(x = deathin, order = c(0, 1, 1), include.mean = FALSE)
##
## Coefficients:
## ma1
## -0.8431
## s.e. 0.0422
##
## sigma^2 estimated as 27217: log likelihood = -770.54, aic = 1545.08
predict(arma,7)
## $pred
## Time Series:
## Start = c(28, 1)
## End = c(28, 7)
## Frequency = 7
## [1] 487.4141 487.4141 487.4141 487.4141 487.4141 487.4141 487.4141
##
## $se
## Time Series:
## Start = c(28, 1)
## End = c(28, 7)
## Frequency = 7
## [1] 164.9756 166.9950 168.9903 170.9623 172.9118 174.8395 176.7463