Time Series Decomposition for Houses Sold in the US

1 Introduction

The following data is from https://www.census.gov. The data set contains the number of homes sold (by the thousands) in the United States from the years 1963 to 2024. Each year contains 12 observations for each month of the year with the exception of 2024 that only has 9 observations. There are a total of 741 observations. For this analysis, we will only use the 175 most recent observations.

1.1 Variable Description

• Period - The year and month the number of houses sold was recorded.
• Value - The number of houses bought in the United states in thousands. (Example: 400 = 400,000 houses sold)

1.2 Practical Question

What patterns can be found in our time series and what methods of smoothing and decomposition should be used on it?

2 Exploratory Analysis

First, the data is downloaded. Next, the data is cut down to the 175 most recent observations. this means we are only using observations from the years 2010 to 2024.We also modify the Vale variable so it does not contain commas. There appears to be no missing values in the data.

house <- read.csv("https://raw.githubusercontent.com/AvaDeSt/STA-321/refs/heads/main/Book%203(Sheet1).csv", header = TRUE)
n.row = dim(house)[1]
data.house = house[(n.row-150):n.row, ]
          
data.house$Value = as.numeric(gsub(",", "", data.house$Value))

2.1 Define Time Series Object

Here is a graph of the data with the number of houses sold by the thousands on the y axis. Since the data is recorded by month. Frequency = 12 will be used to define the time series object.

house.ts = ts(data.house[,2], frequency = 12, start = c(2012, 1))
par(mar=c(2,2,2,2))
plot(house.ts, main="US Houses Bought Between March, 2012 and Sept, 2024", ylab="Houses Bought", xlab="")

US Houses Bought Monthly

We can see that the number of houses sold increases gradually throughout the years with a large spike right after 2020. The time series for the most part is additive.

3 Forecasting with Decomposing

Next, we will try forecasting using both classical decomposing and STL decomposing to see which one performs better.

cls.decomp = decompose(house.ts)
par(mar=c(2,2,2,2))
plot(cls.decomp, xlab="")

Classical decomposition of additive time series

stl.decomp=stl(house.ts, s.window = 12)
par(mar=c(2,2,2,2))
plot(stl.decomp)

STL decomposition of additive time series

We can see that in this case, STL decomposition works a little better for our data. While the time series is additive and relatively constant, there is a large spike in the number of houses sold right after 2020, creating a few outliers. Since STL decomposition is more robust, it works better in this case.

4 Data Training and Testing

Next, we will hold the last ten periods of the data for testing. We will define four different training data sets. The training set sizes used in this analysis are 144, 109, 73, and 48. The same test set with size ten will be used to calculate the prediction error.

ini.data = data.house[,2]
n0 = length(ini.data)
##
train.data01 = data.house[1:(n0-7), 2]
train.data02 = data.house[37:(n0-7), 2]
train.data03 = data.house[73:(n0-7), 2]
train.data04 = data.house[97:(n0-7), 2]
## last 7 observations
test.data = data.house[(n0-6):n0,2]
##
train01.ts = ts(train.data01, frequency = 12, start = c(2012, 1))
train02.ts = ts(train.data02, frequency = 12, start = c(2015, 1))
train03.ts = ts(train.data03, frequency = 12, start = c(2018, 1))
train04.ts = ts(train.data04, frequency = 12, start = c(2020, 1))
##
stl01 = stl(train01.ts, s.window = 12)
stl02 = stl(train02.ts, s.window = 12)
stl03 = stl(train03.ts, s.window = 12)
stl04 = stl(train04.ts, s.window = 12)
## Forecast with decomposing
fcst01 = forecast(stl01,h=10, method="naive")
fcst02 = forecast(stl02,h=10, method="naive")
fcst03 = forecast(stl03,h=10, method="naive")
fcst04 = forecast(stl04,h=10, method="naive")

Next, we perform error analysis.

PE01=(test.data-fcst01$mean)/fcst01$mean

## Warning in `-.default`(test.data, fcst01$mean): longer object length is not a
## multiple of shorter object length

PE02=(test.data-fcst02$mean)/fcst02$mean

## Warning in `-.default`(test.data, fcst02$mean): longer object length is not a
## multiple of shorter object length

PE03=(test.data-fcst03$mean)/fcst03$mean

## Warning in `-.default`(test.data, fcst03$mean): longer object length is not a
## multiple of shorter object length

PE04=(test.data-fcst04$mean)/fcst04$mean

## Warning in `-.default`(test.data, fcst04$mean): longer object length is not a
## multiple of shorter object length

###
MAPE1 = mean(abs(PE01))
MAPE2 = mean(abs(PE02))
MAPE3 = mean(abs(PE03))
MAPE4 = mean(abs(PE04))
###
E1=test.data-fcst01$mean

## Warning in `-.default`(test.data, fcst01$mean): longer object length is not a
## multiple of shorter object length

E2=test.data-fcst02$mean

## Warning in `-.default`(test.data, fcst02$mean): longer object length is not a
## multiple of shorter object length

E3=test.data-fcst03$mean

## Warning in `-.default`(test.data, fcst03$mean): longer object length is not a
## multiple of shorter object length

E4=test.data-fcst04$mean

## Warning in `-.default`(test.data, fcst04$mean): longer object length is not a
## multiple of shorter object length

##
MSE1=mean(E1^2)
MSE2=mean(E2^2)
MSE3=mean(E3^2)
MSE4=mean(E4^2)
###
MSE=c(MSE1, MSE2, MSE3, MSE4)
MAPE=c(MAPE1, MAPE2, MAPE3, MAPE4)
accuracy=cbind(MSE=MSE, MAPE=MAPE)
row.names(accuracy)=c("n.144", "n.109", "n. 73", "n. 48")
kable(accuracy, caption="Error comparison between forecast results with different sample sizes")

Error comparison between forecast results with different sample sizes

	MSE	MAPE
n.144	3528.487	0.0802875
n.109	3491.738	0.0795013
n. 73	3965.483	0.0852145
n. 48	5327.657	0.1013145

We can see from the table above that a training size of 109 performs the best and has the lowest errors. While the mean square errors look normal, the mean absolute percentage error is well into the thousands. Some possible reasons for this are the observations in our time series ranges from about 300 to a little over one thousand. The time series also follows several patterns including seasonal trends and being additive. We will take a closer look at the errors next by making one graph for the MSE and one for the MAPE.

par(mfrow=c(1,2))
plot(1:4, MSE, type="b", col="darkred", ylab="Error", xlab="",
     #ylim=c(0.4,.85),xlim = c(0.5,4.5), 
     main="MSE", axes=FALSE)
labs=c("n=144", "n=109", "n=73", "n=48")
axis(1, at=1:4, label=labs)
axis(2)
#lines(1:4, MAPE, type="b", col="blue")
text(1:4, MAPE+0.03, as.character(round(MAPE,4)), col="blue", cex=0.7)
text(1:4, MSE-0.03, as.character(round(MSE,4)), col="darkred", cex=0.7)
legend(1.5, 0.63, c("MSE", "MAPE"), col=c("darkred","blue"), lty=1, bty="n", cex=0.7)
###
#```{r fig.align='center', fig.cap= "Comparing forecast errors", fig.width=5, fig.height=3.5}
plot(1:4, MAPE, type="b", col="darkred", ylab="Error", xlab="",
     #ylim=c(0.4,.85),xlim = c(0.5,4.5), 
     main="MAPE", axes=FALSE)
labs=c("n=144", "n=109", "n=73", "n=48")
axis(1, at=1:4, label=labs)
axis(2)

Comparing forecast errors

5 Summary and Conclusion

To conclude, we examined a time series from the years 2012 to 2024 totaling 151 observations. The data looks at the number of houses sold in the United States for each month of the years recorded. We noticed that the time series follows a mostly seasonal trend that is additive. The only exception was a spike in the observations after 2020 that eventually went back down. We found that STL decomposition works better for this time series because it is more robust. Were this analysis to be continued in the future, I would focus on reducing the MAPE. A way I would try to do this would be by using a different model that accounts more for seasonal patterns and possibly outliers that do not follow a pattern since this time series contains a few.