library(fma)
## Loading required package: tseries
## Loading required package: forecast
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: timeDate
## This is forecast 7.2
summary(dole)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4742 20400 56760 221200 392600 856500
par(mfrow=c(2,2))
plot(dole, main="Monthly unemployment benefits")
plot(log(dole), main="Log transformed unemployment benefits")
lambda <- BoxCox.lambda(dole)
lambda
## [1] 0.3290922
plot(BoxCox(dole,lambda), main="Box-Cox transformed unemployment benefits")
summary(usdeaths)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6892 8089 8728 8788 9323 11320
par(mfrow=c(2,2))
plot(usdeaths, main="Monthly Accidental Deaths in the US")
plot(log(usdeaths), main="Log Transformed Accidental Deaths")
lambda2 <- BoxCox.lambda(usdeaths)
plot(BoxCox(usdeaths,lambda2), main = "Box-Cox transformed Accidental Deaths")
plot(sqrt(usdeaths), main = "SQRT transformed Accidental Deaths")
lambda2
## [1] -0.03363775
summary(bricksq)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 187.0 338.5 428.0 408.8 490.5 589.0
par(mfrow=c(2,2))
plot(bricksq, main="Quarterly Brick Production")
plot(log(bricksq), main = "Log Transformed Brick Production")
lambda3 <- BoxCox.lambda(bricksq)
lambda3
## [1] 0.2548929
plot(BoxCox(bricksq,lambda3), main = "Box-Cox Brick Production")
plot(sqrt(bricksq), main = "SQRT transformed Brick Production")
***
par(mfrow=c(2,2))
plot(dowjones, main="Dow Jones Index")
plot(log(dowjones), main = "Log DJI")
plot(sqrt(dowjones), main = "SQRT")
dj_drift <- rwf(dowjones , h=10, drift=TRUE)
dj_drift_log <- rwf(log(dowjones), h = 10, drift = TRUE)
dj_drift_sqrt <- rwf(sqrt(dowjones), h =10, drift = TRUE)
par(mfrow=c(2,2))
plot(dj_drift,plot.conf=FALSE,main="Drift Method Dow Jones", ylab="Index",xlab="Year")
legend("topleft",lty=1, col=c(4),legend=c("SQRT"))
plot(dj_drift_log,plot.conf=FALSE,main="Log Method Dow Jones", ylab="Index",xlab="Year")
legend("topleft",lty=1, col=c(4),legend=c("Log"))
plot(dj_drift_sqrt,plot.conf=FALSE,main="SQRT Method Dow Jones", ylab="Index",xlab="Year")
legend("topleft",lty=1, col=c(4),legend=c("Drift"))
dj2 <- window(dowjones, start=1, end=66-.1)
dj2a <- meanf(dj2,h=12)
dj2b <- rwf(dj2,h=12)
dj2c <- rwf(dj2,h=12, drift = TRUE)
plot(dj2a, plot.conf=FALSE, main="Dow Jones Index", xlim=c(1,78))
lines(dj2b$mean,col=2)
lines(dj2c$mean,col=3)
lines(dowjones)
legend("topleft", lty=1, col=c(4,2,3), legend=c("Mean ","Naive","Drifit"))
dj_drift <- rwf(dowjones , h=24, drift=TRUE)
dj_mean <-meanf(dowjones, h=42)
dj_naive <-naive(dowjones, h=42)
plot(dj_drift,plot.conf=FALSE,main="Drift Method Dow Jones", ylab="Index",xlab="Year")
lines(dj_mean$mean, col=2)
lines(dj_naive$mean, col=3)
legend("topleft",lty=1, col=c(4,2,3),legend=c("Mean Method","Naive Method","Drift"))
head(ibmclose)
## [1] 460 457 452 459 462 459
summary(ibmclose)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 306.0 387.0 494.0 478.5 549.0 603.0
par(mfrow=c(2,2))
plot(ibmclose)
qqnorm(ibmclose)
qqline(ibmclose)
plot(log(ibmclose))
plot(sqrt(ibmclose))
ibmclose_train <- window(ibmclose ,end=300)
ibmclose_test <- window(ibmclose ,start=301)
ibm_avg <- meanf(ibmclose_train,h=54)$mean
ibm_naive <- naive(ibmclose_train ,h=54)$mean
ibm_drift <- rwf(ibmclose_train ,drift=TRUE,h=54)$mean
plot(ibmclose_train,main="IBM Close Prices",xlab="Day",ylab="Price")
lines(ibm_naive,col=2)
lines(ibm_avg,col=4)
lines(ibm_drift,col=3)
lines(ibmclose_test,col=8)
legend("topleft",lty=1,col=c(4,2,3),
legend=c("Mean Method","Naive Method","Drift Method"))
plot(ibmclose_train,main="IBM Close Prices", ylab="Price",xlab="Day", xlim=c(250,369), ylim=c(300,425))
lines(ibm_naive,col=2)
lines(ibm_avg,col=4)
lines(ibm_drift,col=3)
lines(ibmclose_test,col=8)
legend("topleft",lty=1,col=c(4,2,3),
legend=c("Mean Method","Naive Method","Drift Method"))
head(hsales)
## [1] 55 60 68 63 65 61
summary(hsales)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 24.00 44.00 53.00 52.29 59.00 89.00
par(mfrow=c(2,2))
plot(hsales)
qqnorm(hsales)
qqline(hsales)
plot(log(hsales))
acf(hsales)
hsales_ts <- ts(hsales,start=1,end=275)
hsales_train <- window(hsales_ts,end=251)
hsales_test <- window(hsales_ts,start=251)
hsales_avg <- meanf(hsales_train,h=24)$mean
hsales_naive <- naive(hsales_train,h=24)$mean
hsales_drift <- rwf(hsales_train,drift=TRUE,h=24)$mean
plot(hsales_train,main="House Sales",xlab="Month",ylab="Price")
lines(hsales_naive,col=2)
lines(hsales_avg,col=4)
lines(hsales_drift,col=3)
lines(hsales_test,col=8)
legend("topleft",lty=1,col=c(4,2,3),
legend=c("Mean Method","Naive Method","Drift Method"))
plot(hsales_train,main="House sales", ylab="House Sales",xlab="Month", xlim=c(240,275), ylim=c(35,75))
lines(hsales_naive,col=2)
lines(hsales_avg,col=4)
lines(hsales_drift,col=3)
lines(hsales_test,col=8)
legend("topleft",lty=1,col=c(4,2,3), legend=c("Mean Method","Naive Method","Drift Method"))
econsumption
## Mwh temp
## 1 16.3 29.3
## 2 16.8 21.7
## 3 15.5 23.7
## 4 18.2 10.4
## 5 15.2 29.7
## 6 17.5 11.9
## 7 19.8 9.0
## 8 19.0 23.4
## 9 17.5 17.8
## 10 16.0 30.0
## 11 19.6 8.6
## 12 18.0 11.8
plot(Mwh ~ temp, data = econsumption, main = "Econsumption")
fit = lm(formula = Mwh ~ temp, data = econsumption)
abline(fit, col=5)
summary(fit)
##
## Call:
## lm(formula = Mwh ~ temp, data = econsumption)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2593 -0.5395 -0.1827 0.4274 2.1972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.19952 0.73040 27.66 8.86e-11 ***
## temp -0.14516 0.03549 -4.09 0.00218 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9888 on 10 degrees of freedom
## Multiple R-squared: 0.6258, Adjusted R-squared: 0.5884
## F-statistic: 16.73 on 1 and 10 DF, p-value: 0.00218
par(mfrow=c(2,2))
plot(fit)
coeffs = coefficients(fit)
pred_temp = c(10, 35)
p_temp = coeffs[1] + coeffs[2]*pred_temp
p_temp
## [1] 18.74795 15.11902
par(mfrow=c(1,2))
fcast <- forecast(fit, newdata=data.frame(temp=10))
plot(fcast, xlab="temp", ylab="Mwh")
fcast2 <- forecast(fit, newdata=data.frame(temp=35))
plot(fcast2, xlab="temp", ylab="Mwh")
temp10 = data.frame(temp=10)
temp35 = data.frame(temp=35)
predict(fit, temp10, interval="predict")
## fit lwr upr
## 1 18.74795 16.34824 21.14766
predict(fit, temp35, interval="predict")
## fit lwr upr
## 1 15.11902 12.49768 17.74035
olympic
).olympic
to include the winning times from the last few Olympics.olympic1 <- matrix(c(1896, 54.2, 1900, 49.4, 1904, 49.2, 1908, 50, 1912 , 48.2, 1920 , 49.6, 1924 , 47.6 , 1928 , 47.8, 1932, 46.2, 1936, 46.5, 1948, 46.2, 1952, 45.9, 1956, 46.7, 1960, 44.9, 1964, 45.1, 1968, 43.8 , 1972, 44.66, 1976, 44.26, 1980, 44.6, 1984, 44.27, 1988 , 43.87 , 1992, 43.5, 1996 , 43.49 , 2000, 43.84 , 2004, 44, 2008, 43.75, 2012, 43.94, 2016 , 43.03) ,ncol=2,byrow=TRUE)
colnames(olympic1) <- c("Year","time")
olympic_ts <- ts(olympic1,start=1,end=28)
plot(time ~ Year, data = olympic_ts, main = "Olympic Gold Medal Times")
fit1 = lm(formula = time ~ Year, data = olympic_ts)
plot(time ~ Year, data = olympic_ts, main = "Olympic Gold Medal Times")
abline(fit1, col=5)
summary(fit1)
##
## Call:
## lm(formula = time ~ Year, data = olympic_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6002 -0.5747 -0.2858 0.5751 4.1505
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 172.481477 11.487522 15.02 2.52e-14 ***
## Year -0.064574 0.005865 -11.01 2.75e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.136 on 26 degrees of freedom
## Multiple R-squared: 0.8234, Adjusted R-squared: 0.8166
## F-statistic: 121.2 on 1 and 26 DF, p-value: 2.752e-11
par(mfrow=c(2,2))
plot(fit1)
coeffs1 = coefficients(fit1)
pred_time = c(2000, 2004, 2008, 2012)
p_time = coeffs1[1] + coeffs1[2]*pred_time
p_time
## [1] 43.33379 43.07549 42.81720 42.55890
par(mfrow=c(2,2))
fcast3 <- forecast(fit1, newdata=data.frame(Year=2000))
plot(fcast3, xlab="Year", ylab="time")
fcast4 <- forecast(fit1, newdata=data.frame(Year=2004))
plot(fcast4, xlab="Year", ylab="time")
fcast5 <- forecast(fit1, newdata=data.frame(Year=2008))
plot(fcast5, xlab="Year", ylab="time")
fcast6 <- forecast(fit1, newdata=data.frame(Year=2012))
plot(fcast6, xlab="Year", ylab="time")
time2000 = data.frame(Year=2000)
time2004 = data.frame(Year=2004)
time2008 = data.frame(Year=2008)
time2012 = data.frame(Year=2012)
predict(fit1, time2000, interval="predict")
## fit lwr upr
## 1 43.33379 40.90529 45.76229
predict(fit1, time2004, interval="predict")
## fit lwr upr
## 1 43.07549 40.63659 45.5144
predict(fit1, time2008, interval="predict")
## fit lwr upr
## 1 42.8172 40.36698 45.26741
predict(fit1, time2012, interval="predict")
## fit lwr upr
## 1 42.5589 40.09648 45.02132
coeffs1 = coefficients(fit1)
pred_time = c(2000, 2004, 2008, 2012)
p_time = coeffs1[1] + coeffs1[2]*pred_time
p_time
## [1] 43.33379 43.07549 42.81720 42.55890
olympic_ts[24:27,2]
## [1] 43.84 44.00 43.75 43.94
logy=B0+B1logx+e
take differential, assuming other variables are constant
dy/y = dx/x*B1 (dy/y / dx/x) = B1
change in y divided by y over the the change in x divided by x
percentage in y over the percentage change in x B1 is elasticity- the partial elasticity of the dependent variable with respect to the independent variable, ie percentage increase in the dependent variable with the percentage increase in the independent variable.
100(dy/y) = 100(dx/x)B1 % change y = % change x B1