Original Model and Preliminary Stuff
Cow <-read.csv("C:/Users/caitl_o8hbexg/Downloads/FinalTS.csv")
MV <-read.csv("C:/Users/caitl_o8hbexg/Downloads/MMM.csv")
attach(Cow)
plot(Milk.Volume~Year, main="Pounds of Milk per Cow")
plot.ts(Milk.Volume, main="Pounds of Milk per Cow")
mod<-lm(Milk.Volume~Year)
summary(mod)
##
## Call:
## lm(formula = Milk.Volume ~ Year)
##
## Residuals:
## Min 1Q Median 3Q Max
## -101.04 -50.02 -15.30 42.88 139.05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -39239.40 2283.24 -17.19 <2e-16 ***
## Year 20.31 1.16 17.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 60.74 on 166 degrees of freedom
## Multiple R-squared: 0.6489, Adjusted R-squared: 0.6468
## F-statistic: 306.8 on 1 and 166 DF, p-value: < 2.2e-16
mod1<-lm(Milk.Volume~Year+I(Year^2))
summary(mod1)
##
## Call:
## lm(formula = Milk.Volume ~ Year + I(Year^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -92.95 -51.66 -15.73 44.30 134.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.234e+06 1.236e+06 -1.808 0.0724 .
## Year 2.250e+03 1.255e+03 1.792 0.0749 .
## I(Year^2) -5.661e-01 3.187e-01 -1.776 0.0776 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 60.35 on 165 degrees of freedom
## Multiple R-squared: 0.6555, Adjusted R-squared: 0.6513
## F-statistic: 157 on 2 and 165 DF, p-value: < 2.2e-16
plot.ts(Milk.Volume, main = "Quadratic Model")
abline(-39237.72, 20.31)
forecast<-coef(summary(mod))[1,1]+coef(summary(mod))[2,1]*1980
forecast
## [1] 977.2886
par(mfrow=c(2,1))
plot(Milk.Volume~Year)
with(Cow, plot(Milk.Volume~Year, type="l"))
Adding Month to Model
justyear<-floor(Cow$Year)
m <- Cow$Year - justyear
mfactor <-factor(round(m*12))
cbind(Cow$Year, mfactor)
## mfactor
## [1,] 1962.083 2
## [2,] 1962.167 3
## [3,] 1962.250 4
## [4,] 1962.333 5
## [5,] 1962.417 6
## [6,] 1962.500 7
## [7,] 1962.583 8
## [8,] 1962.667 9
## [9,] 1962.750 10
## [10,] 1962.833 11
## [11,] 1962.917 12
## [12,] 1963.000 13
## [13,] 1963.083 2
## [14,] 1963.167 3
## [15,] 1963.250 4
## [16,] 1963.333 5
## [17,] 1963.417 6
## [18,] 1963.500 7
## [19,] 1963.583 8
## [20,] 1963.667 9
## [21,] 1963.750 10
## [22,] 1963.833 11
## [23,] 1963.917 12
## [24,] 1964.000 13
## [25,] 1964.083 2
## [26,] 1964.167 3
## [27,] 1964.250 4
## [28,] 1964.333 5
## [29,] 1964.417 6
## [30,] 1964.500 7
## [31,] 1964.583 8
## [32,] 1964.667 9
## [33,] 1964.750 10
## [34,] 1964.833 11
## [35,] 1964.917 12
## [36,] 1965.000 13
## [37,] 1965.083 2
## [38,] 1965.167 3
## [39,] 1965.250 4
## [40,] 1965.333 5
## [41,] 1965.417 6
## [42,] 1965.500 7
## [43,] 1965.583 8
## [44,] 1965.667 9
## [45,] 1965.750 10
## [46,] 1965.833 11
## [47,] 1965.917 12
## [48,] 1966.000 13
## [49,] 1966.083 2
## [50,] 1966.167 3
## [51,] 1966.250 4
## [52,] 1966.333 5
## [53,] 1966.417 6
## [54,] 1966.500 7
## [55,] 1966.583 8
## [56,] 1966.667 9
## [57,] 1966.750 10
## [58,] 1966.833 11
## [59,] 1966.917 12
## [60,] 1967.000 13
## [61,] 1967.083 2
## [62,] 1967.167 3
## [63,] 1967.250 4
## [64,] 1967.333 5
## [65,] 1967.417 6
## [66,] 1967.500 7
## [67,] 1967.583 8
## [68,] 1967.667 9
## [69,] 1967.750 10
## [70,] 1967.833 11
## [71,] 1967.917 12
## [72,] 1968.000 13
## [73,] 1968.083 2
## [74,] 1968.167 3
## [75,] 1968.250 4
## [76,] 1968.333 5
## [77,] 1968.417 6
## [78,] 1968.500 7
## [79,] 1968.583 8
## [80,] 1968.667 9
## [81,] 1968.750 10
## [82,] 1968.833 11
## [83,] 1968.917 12
## [84,] 1969.000 13
## [85,] 1969.083 2
## [86,] 1969.167 3
## [87,] 1969.250 4
## [88,] 1969.333 5
## [89,] 1969.417 6
## [90,] 1969.500 7
## [91,] 1969.583 8
## [92,] 1969.667 9
## [93,] 1969.750 10
## [94,] 1969.833 11
## [95,] 1969.917 12
## [96,] 1970.000 13
## [97,] 1970.083 2
## [98,] 1970.167 3
## [99,] 1970.250 4
## [100,] 1970.333 5
## [101,] 1970.417 6
## [102,] 1970.500 7
## [103,] 1970.583 8
## [104,] 1970.667 9
## [105,] 1970.750 10
## [106,] 1970.833 11
## [107,] 1970.917 12
## [108,] 1971.000 13
## [109,] 1971.083 2
## [110,] 1971.167 3
## [111,] 1971.250 4
## [112,] 1971.333 5
## [113,] 1971.417 6
## [114,] 1971.500 7
## [115,] 1971.583 8
## [116,] 1971.667 9
## [117,] 1971.750 10
## [118,] 1971.833 11
## [119,] 1971.917 12
## [120,] 1972.000 13
## [121,] 1972.083 2
## [122,] 1972.167 3
## [123,] 1972.250 4
## [124,] 1972.333 5
## [125,] 1972.417 6
## [126,] 1972.500 7
## [127,] 1972.583 8
## [128,] 1972.667 9
## [129,] 1972.750 10
## [130,] 1972.833 11
## [131,] 1972.917 12
## [132,] 1973.000 13
## [133,] 1973.083 2
## [134,] 1973.167 3
## [135,] 1973.250 4
## [136,] 1973.333 5
## [137,] 1973.417 6
## [138,] 1973.500 7
## [139,] 1973.583 8
## [140,] 1973.667 9
## [141,] 1973.750 10
## [142,] 1973.833 11
## [143,] 1973.917 12
## [144,] 1974.000 13
## [145,] 1974.083 2
## [146,] 1974.167 3
## [147,] 1974.250 4
## [148,] 1974.333 5
## [149,] 1974.417 6
## [150,] 1974.500 7
## [151,] 1974.583 8
## [152,] 1974.667 9
## [153,] 1974.750 10
## [154,] 1974.833 11
## [155,] 1974.917 12
## [156,] 1975.000 13
## [157,] 1975.083 2
## [158,] 1975.167 3
## [159,] 1975.250 4
## [160,] 1975.333 5
## [161,] 1975.417 6
## [162,] 1975.500 7
## [163,] 1975.583 8
## [164,] 1975.667 9
## [165,] 1975.750 10
## [166,] 1975.833 11
## [167,] 1975.917 12
## [168,] 1976.000 1
model = lm(Milk.Volume ~ justyear + mfactor)
summary(model)
##
## Call:
## lm(formula = Milk.Volume ~ justyear + mfactor)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.268 -9.116 -0.710 6.297 37.032
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.019e+04 6.165e+02 -65.194 < 2e-16 ***
## justyear 2.077e+01 3.119e-01 66.584 < 2e-16 ***
## mfactor1 3.982e+01 1.689e+01 2.358 0.019646 *
## mfactor2 2.034e+00 1.689e+01 0.120 0.904299
## mfactor3 9.625e+01 1.689e+01 5.699 5.97e-08 ***
## mfactor4 1.128e+02 1.689e+01 6.680 4.11e-10 ***
## mfactor5 1.755e+02 1.689e+01 10.393 < 2e-16 ***
## mfactor6 1.489e+02 1.689e+01 8.816 2.38e-15 ***
## mfactor7 1.008e+02 1.689e+01 5.969 1.58e-08 ***
## mfactor8 6.025e+01 1.689e+01 3.567 0.000481 ***
## mfactor9 1.939e+01 1.689e+01 1.148 0.252694
## mfactor10 2.461e+01 1.689e+01 1.457 0.147189
## mfactor11 -4.680e+00 1.689e+01 -0.277 0.782063
## mfactor12 3.475e+01 1.695e+01 2.049 0.042111 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.16 on 154 degrees of freedom
## Multiple R-squared: 0.9769, Adjusted R-squared: 0.975
## F-statistic: 502 on 13 and 154 DF, p-value: < 2.2e-16
par(mfrow=c(1,1))
plot(Milk.Volume~Year, type = "l")
lines(Cow$Year, model$fitted.values, type = "l", col = "red")
Smoothing
out <- filter(Milk.Volume, sides=2, filter=rep(1/3,3)) # moving average
par(mfrow=c(2,1))
plot.ts(Milk.Volume, main="Unsmoothed")
plot.ts(out, main="Smoothed")
CowTS <- ts(MV, frequency =12 , start = c(1962,1))
k <- ksmooth(time(CowTS), CowTS, kernel = "normal", bandwidth = 2)
par(mfrow=c(1,1))
plot(CowTS)
lines(k, lwd=2, col=6)
low <- lowess(CowTS, f=.04)
plot(CowTS)
lines(low, lwd=2, col=5)
spline <- smooth.spline(time(CowTS), CowTS, spar = .2)
plot(CowTS, main = "Smoothed Spline")
lines(spline, lwd=2, col = 7)
ARIMA
library(quantmod)
## Warning: package 'quantmod' was built under R version 3.4.4
## Loading required package: xts
## Warning: package 'xts' was built under R version 3.4.4
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.4.4
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(tseries)
## Warning: package 'tseries' was built under R version 3.4.4
library(timeSeries)
## Warning: package 'timeSeries' was built under R version 3.4.4
## Loading required package: timeDate
## Warning: package 'timeDate' was built under R version 3.4.4
##
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
##
## time<-
library(forecast)
## Warning: package 'forecast' was built under R version 3.4.4
library(xts)
modA <- auto.arima(CowTS)
modA
## Series: CowTS
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.2204 -0.6214
## s.e. 0.0748 0.0627
##
## sigma^2 estimated as 53.42: log likelihood=-530.15
## AIC=1066.3 AICc=1066.46 BIC=1075.43
fore <- forecast(modA, h=12)
fore
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1976 864.9773 855.6103 874.3443 850.6517 879.3029
## Feb 1976 817.7493 805.8719 829.6267 799.5843 835.9142
## Mar 1976 924.4056 910.4626 938.3485 903.0817 945.7295
## Apr 1976 937.4836 921.7439 953.2233 913.4118 961.5554
## May 1976 1000.6235 983.2721 1017.9749 974.0868 1027.1601
## Jun 1976 973.2165 954.3909 992.0420 944.4252 1002.0077
## Jul 1976 931.8501 911.6576 952.0426 900.9684 962.7318
## Aug 1976 892.2597 870.7873 913.7322 859.4204 925.0991
## Sep 1976 846.3679 823.6875 869.0483 811.6812 881.0545
## Oct 1976 851.5326 827.7055 875.3597 815.0921 887.9731
## Nov 1976 817.4931 792.5719 842.4143 779.3795 855.6068
## Dec 1976 859.7534 833.7842 885.7225 820.0370 899.4698
plot(fore)
ARIMA Diagnostics
TScomp <- decompose(CowTS)
plot(TScomp)
diffseas <- diff(CowTS, lag=12)
plot(diffseas)
adf.test(diffseas)
##
## Augmented Dickey-Fuller Test
##
## data: diffseas
## Dickey-Fuller = -2.8615, Lag order = 5, p-value = 0.2172
## alternative hypothesis: stationary
Acf(diffseas)
Pacf(diffseas)
modB <- arima(CowTS, order = c(1, 1, 6), season = list(order = c(0, 1, 1), period= 12))
summary(modB)
##
## Call:
## arima(x = CowTS, order = c(1, 1, 6), seasonal = list(order = c(0, 1, 1), period = 12))
##
## Coefficients:
## ar1 ma1 ma2 ma3 ma4 ma5 ma6 sma1
## -0.8059 0.5835 -0.1454 0.1644 -0.027 -0.1486 -0.1225 -0.6089
## s.e. 0.2802 0.2916 0.1239 0.0988 0.103 0.1123 0.1205 0.0705
##
## sigma^2 estimated as 50.17: log likelihood = -526.65, aic = 1071.31
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01048012 6.805862 4.931041 -0.005252122 0.6582728 0.1264175
## ACF1
## Training set 0.005888113
teststats <- coef(modB)/sqrt(diag(vcov(modB)))
2*pt(abs(teststats), df=160, lower.tail=FALSE)
## ar1 ma1 ma2 ma3 ma4
## 4.567448e-03 4.707633e-02 2.425208e-01 9.796569e-02 7.932051e-01
## ma5 ma6 sma1
## 1.877185e-01 3.108260e-01 5.682440e-15
modC <- arima(CowTS, order = c(1, 1, 5), season = list(order = c(0, 1, 1), period= 12))
teststats <- coef(modC)/sqrt(diag(vcov(modC)))
2*pt(abs(teststats), df=161, lower.tail=FALSE)
## ar1 ma1 ma2 ma3 ma4
## 5.741302e-01 3.482158e-01 4.985974e-01 2.322464e-01 1.824110e-01
## ma5 sma1
## 9.084421e-01 1.967194e-17
modD <- arima(CowTS, order = c(1, 1, 4), season = list(order = c(0, 1, 1), period= 12))
teststats <- coef(modD)/sqrt(diag(vcov(modD)))
2*pt(abs(teststats), df=162, lower.tail=FALSE)
## ar1 ma1 ma2 ma3 ma4
## 2.945397e-01 9.790050e-02 3.413833e-01 2.412106e-01 3.595116e-02
## sma1
## 1.623795e-17
modE <- arima(CowTS, order = c(1, 1, 3), season = list(order = c(0, 1, 1), period= 12))
teststats <- coef(modE)/sqrt(diag(vcov(modE)))
2*pt(abs(teststats), df=163, lower.tail=FALSE)
## ar1 ma1 ma2 ma3 sma1
## 1.743726e-12 5.114379e-06 1.676342e-01 1.598048e-01 9.372782e-17
modF <- arima(CowTS, order = c(1, 1, 2), season = list(order = c(0, 1, 1), period= 12))
teststats <- coef(modF)/sqrt(diag(vcov(modF)))
2*pt(abs(teststats), df=164, lower.tail=FALSE)
## ar1 ma1 ma2 sma1
## 6.494291e-01 3.180498e-01 3.472382e-01 3.288593e-18
modG <- arima(CowTS, order = c(1, 1, 1), season = list(order = c(0, 1, 1), period= 12))
teststats <- coef(modG)/sqrt(diag(vcov(modG)))
2*pt(abs(teststats), df=165, lower.tail=FALSE)
## ar1 ma1 sma1
## 6.100689e-01 7.076679e-01 2.072820e-18
modH <- arima(CowTS, order = c(2, 1, 2), season = list(order = c(0, 1, 1), period= 12))
teststats <- coef(modH)/sqrt(diag(vcov(modH)))
2*pt(abs(teststats), df=163, lower.tail=FALSE)
## ar1 ar2 ma1 ma2 sma1
## 1.356270e-01 1.115922e-03 2.854783e-03 7.894490e-06 9.018423e-18
Box.test(modH$residuals, type = "Ljung")
##
## Box-Ljung test
##
## data: modH$residuals
## X-squared = 0.18569, df = 1, p-value = 0.6665
Log Likelihood
modH <- arima(CowTS, order = c(2, 1, 2), season = list(order = c(0, 1, 1), period= 12))
modA <- auto.arima(CowTS)
modA
## Series: CowTS
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.2204 -0.6214
## s.e. 0.0748 0.0627
##
## sigma^2 estimated as 53.42: log likelihood=-530.15
## AIC=1066.3 AICc=1066.46 BIC=1075.43
AIC(modH)
## [1] 1068.471
AIC(modA)
## [1] 1066.301
modJ <- arima(CowTS, order = c(0, 1, 1), season = list(order = c(0, 1, 0), period= 12))
modJ
##
## Call:
## arima(x = CowTS, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 0), period = 12))
##
## Coefficients:
## ma1
## -0.2052
## s.e. 0.0749
##
## sigma^2 estimated as 75.91: log likelihood = -555.5, aic = 1114.99
logLik(modA)
## 'log Lik.' -530.1504 (df=3)
logLik(modJ)
## 'log Lik.' -555.497 (df=2)
t1<-as.numeric(logLik(modA))
t2<-as.numeric(logLik(modJ))
teststat<- 2*(t1-t2)
pchisq(teststat, df=1, lower.tail = FALSE)
## [1] 1.079946e-12
finalforecast <- forecast(modA, h=12)
finalforecast
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1976 864.9773 855.6103 874.3443 850.6517 879.3029
## Feb 1976 817.7493 805.8719 829.6267 799.5843 835.9142
## Mar 1976 924.4056 910.4626 938.3485 903.0817 945.7295
## Apr 1976 937.4836 921.7439 953.2233 913.4118 961.5554
## May 1976 1000.6235 983.2721 1017.9749 974.0868 1027.1601
## Jun 1976 973.2165 954.3909 992.0420 944.4252 1002.0077
## Jul 1976 931.8501 911.6576 952.0426 900.9684 962.7318
## Aug 1976 892.2597 870.7873 913.7322 859.4204 925.0991
## Sep 1976 846.3679 823.6875 869.0483 811.6812 881.0545
## Oct 1976 851.5326 827.7055 875.3597 815.0921 887.9731
## Nov 1976 817.4931 792.5719 842.4143 779.3795 855.6068
## Dec 1976 859.7534 833.7842 885.7225 820.0370 899.4698
plot(finalforecast)