Stat lab 15
Rebecca Lewis
April 10, 2018
STAT 360 Lab 15: Forecasting
Name: Rebecca Lewis
Set Up Work Space
Load R libraries
library(rmarkdown)
library(knitr)
library(forecast)Set the Seed
set.seed(33)Load the Kansas City property values data
propertyValues <- read.csv(url("https://www.dropbox.com/s/kbzr00qy4b9kks3/STAT_360_-_Property_Values.csv?dl=1"))
attach(propertyValues)Visualizing a Time Series
Exercise 1:
Use tapply() to calculate the mean Selling Price of properties as a function of Year Built (i.e., the Selling Price of properties built in 1900, the mean Selling Price of properties built in 1901, etc.), assign it to the object averagedPrices, and display the resulting values. Hint: the tapply() function takes three arguments - a vector of the original values, a vector of delineations by which to separate the original values, and the name of the function that should be applied to each delineated group (e.g., mean, median, max).
averagedPrices<- tapply(SellingPrice, YearBuilt, mean)
averagedPrices## 1900 1901 1902 1903 1904 1905 1906 1907
## 581536.6 557108.3 673192.6 480958.2 583867.8 753443.9 670027.7 676324.5
## 1908 1909 1910 1911 1912 1913 1914 1915
## 564499.8 696449.0 671671.8 632584.2 613193.2 586066.3 615246.1 585036.9
## 1916 1917 1918 1919 1920 1921 1922 1923
## 601041.6 528126.8 492346.9 537887.6 477761.0 613224.2 569794.1 618653.8
## 1924 1925 1926 1927 1928 1929 1930 1931
## 570419.9 607316.6 625443.4 654154.2 621920.2 574396.8 600659.3 661781.1
## 1932 1933 1934 1935 1936 1937 1938 1939
## 458488.4 772483.3 517152.4 542229.8 642335.0 648791.5 548884.9 585718.4
## 1940 1941 1942 1943 1944 1945 1946 1947
## 576255.1 526657.0 379312.9 333347.0 356322.1 447466.2 524712.0 450121.2
## 1948 1949 1950 1951 1952 1953 1954 1955
## 433355.6 473352.6 490592.9 545169.0 530501.0 490505.5 453709.5 450762.0
## 1956 1957 1958 1959 1960 1961 1962 1963
## 470342.4 484821.1 474179.5 449222.6 453246.6 432230.9 436771.4 493211.7
## 1964 1965 1966 1967 1968 1969 1970 1971
## 514325.8 496542.7 446796.4 444060.9 458030.9 425526.5 439935.4 441974.6
## 1972 1973 1974 1975 1976 1977 1978 1979
## 528653.3 541459.7 503944.4 507325.5 526411.7 495485.2 472326.5 486019.1
## 1980 1981 1982 1983 1984 1985 1986 1987
## 491796.0 471219.5 562069.0 505468.3 555777.9 511859.7 476989.1 517565.0
## 1988 1989 1990 1991 1992 1993 1994 1995
## 583930.4 583063.4 564133.4 630630.6 548205.9 556760.5 486864.0 577933.8
## 1996 1997 1998 1999 2000 2001 2002 2003
## 639673.5 606173.9 594280.1 640431.2 682003.6 741340.0 578818.5 558791.4
## 2004 2005 2006 2007 2008 2009 2010 2011
## 596095.0 580895.5 631041.5 615193.3 642037.7 518462.2 551678.4 544648.4
## 2012 2013 2014 2015
## 527437.0 678599.6 683792.7 759970.9Exercise 2:
Use plot() generate a scatter plot of averagedPrices against YearBuilt where the data points are displayed via connected lines. Hint: the plot() argument ‘type’ will come in handy.
year<-c(1900:2015)
plot(averagedPrices~year,type="l", data= propertyValues, main="Average Prices of Homes", xlab="Year", ylab="Average Price of Home")
Exercise 3:
Use the scatter plot to describe the association between the average Selling Price and Year Built.
The average selling price of houses does not seem to be affected by the Year the house was built. There was a massive decrease in price through in the year 1949.Exercise 4:
Use ma() to calculate a fifth-order moving average of averagedPrices, assign it to the object movingPrices, and display the resulting values.
movingPrices<-ma(averagedPrices, order=5)
movingPrices## Time Series:
## Start = 1
## End = 116
## Frequency = 1
## [1] NA NA 575332.7 609714.2 632298.0 632924.4 649632.7
## [8] 672149.0 655794.6 648305.9 635679.6 639992.9 623752.3 606425.3
## [15] 600116.8 583103.5 564359.7 548888.0 527432.8 529869.3 538202.8
## [22] 563464.1 569970.6 595881.7 598325.6 615197.6 615850.9 616646.2
## [29] 615314.8 622582.3 583449.2 613561.8 602112.9 590427.0 586537.8
## [36] 624598.4 579878.7 593591.9 600397.0 577261.4 523365.7 480258.1
## [43] 434378.8 408621.1 408232.0 422393.7 442395.4 465801.5 474426.9
## [50] 478518.3 494594.2 506024.2 502095.6 494129.4 479164.1 470028.1
## [57] 466762.9 465865.5 466362.4 458740.1 449130.2 452936.6 465957.3
## [64] 474616.5 477529.6 478987.5 471951.3 454191.5 442870.0 441905.6
## [71] 458824.1 475509.9 491193.5 504671.5 521558.9 514925.3 501098.7
## [78] 497513.6 494407.7 483369.3 496686.0 503314.4 517266.1 521278.9
## [85] 522432.8 513532.0 529224.4 534681.5 545136.3 575864.6 581992.8
## [92] 576558.8 557318.9 560079.0 561887.5 573481.1 580985.1 611698.5
## [99] 632512.5 652845.8 647374.7 640276.9 631409.7 611188.1 589128.4
## [106] 596403.3 613052.6 597526.0 591682.6 574404.0 556852.7 564165.1
## [113] 597231.2 638889.7 NA NAExercise 5:
Use lines() to overlay movingPrices on top of the scatter plot you generated in Exericse 2 using a thick line. Hint: you will need to re-generate the scatter plot before overlaying movingPrices on top of it.
g_range1 <- range(0,averagedPrices )
plot(averagedPrices,type="l", main="Average Prices of Homes", xlab="Year", ylab="Average Price of Home",axes=FALSE)
axis(1, at=1:117, lab=c(1900:2016))
axis(2, las=1, at=50000*0:g_range1[2])
lines(movingPrices,lwd=4)
Exercise 6:
Do you observe any cyclical seasonality in the association?
There was clear cyclical seasonality in the association, they are moving up and down around the same times.Constructing an ARIMA Model
Exercise 7:
Use ts() to convert the first 101 averagedPrices into a time series spanning the years between 1900 and 2000, assign it to the object timeSeries, and display the result.
timesSeries<-ts(averagedPrices[1:101])
timesSeries## Time Series:
## Start = 1
## End = 101
## Frequency = 1
## 1900 1901 1902 1903 1904 1905 1906 1907
## 581536.6 557108.3 673192.6 480958.2 583867.8 753443.9 670027.7 676324.5
## 1908 1909 1910 1911 1912 1913 1914 1915
## 564499.8 696449.0 671671.8 632584.2 613193.2 586066.3 615246.1 585036.9
## 1916 1917 1918 1919 1920 1921 1922 1923
## 601041.6 528126.8 492346.9 537887.6 477761.0 613224.2 569794.1 618653.8
## 1924 1925 1926 1927 1928 1929 1930 1931
## 570419.9 607316.6 625443.4 654154.2 621920.2 574396.8 600659.3 661781.1
## 1932 1933 1934 1935 1936 1937 1938 1939
## 458488.4 772483.3 517152.4 542229.8 642335.0 648791.5 548884.9 585718.4
## 1940 1941 1942 1943 1944 1945 1946 1947
## 576255.1 526657.0 379312.9 333347.0 356322.1 447466.2 524712.0 450121.2
## 1948 1949 1950 1951 1952 1953 1954 1955
## 433355.6 473352.6 490592.9 545169.0 530501.0 490505.5 453709.5 450762.0
## 1956 1957 1958 1959 1960 1961 1962 1963
## 470342.4 484821.1 474179.5 449222.6 453246.6 432230.9 436771.4 493211.7
## 1964 1965 1966 1967 1968 1969 1970 1971
## 514325.8 496542.7 446796.4 444060.9 458030.9 425526.5 439935.4 441974.6
## 1972 1973 1974 1975 1976 1977 1978 1979
## 528653.3 541459.7 503944.4 507325.5 526411.7 495485.2 472326.5 486019.1
## 1980 1981 1982 1983 1984 1985 1986 1987
## 491796.0 471219.5 562069.0 505468.3 555777.9 511859.7 476989.1 517565.0
## 1988 1989 1990 1991 1992 1993 1994 1995
## 583930.4 583063.4 564133.4 630630.6 548205.9 556760.5 486864.0 577933.8
## 1996 1997 1998 1999 2000
## 639673.5 606173.9 594280.1 640431.2 682003.6Exercise 8:
Use arima() to construct a (10, 1, 7) ARIMA model of the timeSeries, assign it to arima.10.01.07, and use summary() to display the resulting model. Hint: the arima() argument ‘order’ will come in handy.
arima.10.01.07<-arima(timesSeries,order=c(10,1,7))
summary(arima.10.01.07)##
## Call:
## arima(x = timesSeries, order = c(10, 1, 7))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.2300 -0.0811 -0.4130 -0.0444 -0.1326 -0.6598 0.3104 0.2837
## s.e. 0.2301 0.2346 0.2661 0.2706 0.2140 0.2703 0.3042 0.2256
## ar9 ar10 ma1 ma2 ma3 ma4 ma5 ma6
## 0.0781 0.1493 -0.3233 -0.0769 0.3427 -0.2154 0.1869 0.4970
## s.e. 0.2112 0.1589 0.2104 0.1497 0.1701 0.1619 0.1120 0.1725
## ma7
## -0.8890
## s.e. 0.1734
##
## sigma^2 estimated as 3.029e+09: log likelihood = -1239.01, aic = 2514.02
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 672.3583 54762.14 39359.67 -0.8358563 7.410698 0.767515
## ACF1
## Training set 0.001037327Exercise 9:
Use forecast() to predict the next 16 average Selling Prices based on the results of arima.10.01.07, assign them to the object arimaPrediction, and display the results.
arimaPrediction<-forecast(arima.10.01.07,h=16)
arimaPrediction## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 102 649136.2 576850.7 721421.7 538585.1 759687.3
## 103 584813.7 505862.1 663765.2 464067.7 705559.6
## 104 610071.1 525748.4 694393.7 481110.7 739031.4
## 105 598610.6 509914.2 687307.1 462961.1 734260.1
## 106 618372.7 525889.0 710856.4 476931.1 759814.3
## 107 598169.4 500089.3 696249.5 448168.8 748170.0
## 108 613435.7 512520.5 714350.8 459099.2 767772.1
## 109 646425.2 544529.9 748320.6 490589.8 802260.7
## 110 610767.1 506056.7 715477.6 450626.3 770907.9
## 111 609055.0 502255.3 715854.7 445718.9 772391.0
## 112 581357.3 469952.5 692762.2 410978.3 751736.3
## 113 607689.2 495213.8 720164.6 435672.9 779705.5
## 114 593936.1 478147.9 709724.2 416853.4 771018.7
## 115 588276.7 469723.3 706830.1 406964.9 769588.5
## 116 620745.7 500567.1 740924.4 436948.3 804543.1
## 117 619518.5 497955.8 741081.2 433604.4 805432.6Exercise 10:
How does your predicted average Selling Price for 2004 compare to the actual average Selling Price in 2004?
arimaPrediction$mean[4]## [1] 598610.6arimaPrediction$mean[4]-averagedPrices[105]## 2004
## 2515.624The predicted value of the average selling price of a home in 2004 was $598610.6 and the actual value was $596095. This is very close with the prediction being only $2515.624 more than the actual value.Exercise 11:
Use plot() to generate a scatter plot of arimaPrediction and lines() to overlay averagedPrices on top of it.
ar<-arimaPrediction$mean
g_range <- range(0, ar)
plot(ar,type="l", col="blue", axes=FALSE, main="Predicted Average Prices of Homes", xlab="Year", ylab=" Predicted Average Price of Home")
axis(1, at=102:117, lab=c(2001:2016))
axis(2, las=1, at=50000*0:g_range[2])
lines(averagedPrices,lwd=4)
Exercise 12:
How well does arimaPrediction match what actually occured over the next 16 years?
It started out fine predicting the decrease in house prices in the beginning of the 2000s, but past that the prediction did not match the real values since it could not predict the housing price fall in 2009 or the massive price increase in 2013. Exercise 13:
Use arima() to construct a (25, 1, 7) ARIMA model of the timeSeries, use forecast() to predict the next 16 average Selling Prices based on the results of this new model, and finally use plot() to generate a scatter plot of these new predictions and lines() to overlay averagedPrices on top of it
ad<-arima(timesSeries, order=c(25,1,7))
cd<-forecast(ad, h=16)
cd## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 102 641198.6 575389.7 707007.4 540552.6 741844.6
## 103 603198.8 531957.3 674440.3 494244.4 712153.3
## 104 620978.7 544007.7 697949.7 503261.8 738695.6
## 105 577355.7 495277.1 659434.2 451827.4 702884.0
## 106 609281.6 523846.8 694716.4 478620.4 739942.8
## 107 620298.4 529758.6 710838.1 481829.8 758767.0
## 108 633428.3 536358.2 730498.4 484972.4 781884.2
## 109 646063.3 547103.9 745022.6 494718.0 797408.5
## 110 599713.8 499100.5 700327.1 445839.0 753588.6
## 111 608755.9 507021.5 710490.3 453166.6 764345.3
## 112 595638.7 490394.9 700882.5 434682.2 756595.2
## 113 620254.1 514429.0 726079.2 458408.6 782099.7
## 114 610328.9 502258.2 718399.6 445049.0 775608.8
## 115 627698.7 518217.3 737180.1 460261.4 795136.0
## 116 672840.3 561392.8 784287.8 502396.1 843284.5
## 117 636042.7 522463.3 749622.0 462338.0 809747.3ds<-cd$mean
g_range <- range(0, ds)
plot(ds,type="l", col="blue", axes=FALSE,main="Predicted Average Prices of Homes", xlab="Year", ylab=" Predicted Average Price of Home")
axis(1, at=102:117, lab=c(2001:2016))
axis(2, las=1, at=50000*0:g_range[2])
lines(averagedPrices,lwd=4)
Exercise 14:
Use AIC() to determine if the (10, 1, 7) or the (25, 1, 7) ARIMA model of the timeSeries provided a more parsimonious fit.
AIC(ad)## [1] 2530.919AIC(arima.10.01.07)## [1] 2514.023The (10,1,7) of the ARIMA model is a more parsimonious fit.