HousingMarketAnalysis_TS_DL.Rmd

---
title: "Housing Market Analysis - Time Series & Deep Learning"
author: "Sam Vuong, Raymond (shanhua) Huang, Carmon Ho, Kyle Murphy"
date: "July 2020"
output:
  html_document: default
  pdf_document: default
  word_document: default
---

### Abstract 
The goal of this project is to analyze the trend of Melbourne Australia’s housing market. With the housing market on the rise, it began to see a trend of “cool off”. The success of this project is to determine the cause/s of the housing market that has “cool off”. 

Our application on ShinyApps can be accessed at: https://skvuong.shinyapps.io/project/
Our application code can be found in GitHub at: https://github.com/skvuong/housingMarketAnalysis

## Introduction

Every family needs to have a home to be living in, and not everyone can afford to buy a home. There are a lot of variables for a family to be able to secure the home they want and one of the major factors is price. With this project we will examine  Melbourne’s Housing Market on determining factors on why prices have seen decline in January 2016. 
We will employ time series technique to analyze when the housing market has declined.  We will also use the Neural Network Deep Learning technique to model and predict the house price based on other factors.

### Dataset
We are using the “Melbourne Housing Market” from Kaggle website.

The dataset has 34,857 rows and 21 columns. 

The dataset can be found at: https://www.kaggle.com/anthonypino/melbourne-housing-market.

## Install and load the required packages

The following R libraries are used by our application:

```{r, warning=FALSE, error=FALSE, message=FALSE}

# Install tensorflow / keras packages
# Note those will take at least half hours
#install.packages('devtools')
library(devtools)

devtools::install_github("rstudio/tensorflow")
devtools::install_github("rstudio/keras")
#tensorflow::install_tensorflow()
#tensorflow::tf_config()

if(!require(astsa)) 
  install.packages("astsa")
if(!require(dplyr)) 
  install.packages("dplyr")
if(!require(fastDummies)) 
  install.packages("fastDummies")
if(!require(forecast)) 
  install.packages("forecast")
if(!require(corrplot)) 
  install.packages("corrplot")
if(!require(ggplot2)) 
  install.packages("ggplot2")
if(!require(ggthemes)) 
  install.packages("ggthemes")
if(!require(keras)) 
  install.packages("keras")
if(!require(magrittr)) 
  install.packages("magrittr")
if(!require(mlbench)) 
  install.packages("mlbench")
if(!require(neuralnet)) 
  install.packages("neuralnet")
if(!require(RColorBrewer)) 
  install.packages("RColorBrewer")
if(!require(sarima))
  install.packages("sarima")
if(!require(tidyr)) 
  install.packages("tidyr")
if(!require(tseries))
  install.packages("tseries")
if(!require(xts))
  install.packages("xts")

library(astsa)
library(dplyr)
library(fastDummies)
library(forecast)
library(corrplot)
library(ggplot2)
library(ggthemes)
library(keras)
library(magrittr)
library(mlbench)
library(neuralnet)
library(RColorBrewer)
library(sarima)
library(tidyr)
library(tseries)
library(xts)
```

## Loading Data

```{r, warning=FALSE, error=FALSE, message=FALSE}
housing_data <- read.csv("Melbourne_housing_FULL.csv", header=TRUE)
```

## Data Exploration and Visualization

```{r, warning=FALSE, error=FALSE, message=FALSE}
dim(housing_data)

colnames(housing_data)

str(housing_data)

#Number of Unique Suburbs
length(unique(housing_data$Suburb))

#Number of Unique Regions
length(unique(housing_data$Regionname))

#List of Home Types
unique(housing_data$Type)

#1. Top 10 Suburbs By Count of Homes
#Calculate count of homes per suburb
top10 <- housing_data %>% group_by(Suburb) %>%
  summarise(Number = n()) %>% arrange(desc(Number)) %>%
  head(10)
#Plot
ggplot(top10, aes(reorder(Suburb, Number), Number, fill = Suburb))+
  geom_bar(stat = "identity")+
  theme(legend.position = "none")+
  labs(x = "Suburb", y = "Count of Homes",
       title = "Top 10 Suburbs by Count of Homes")+
  coord_flip()

#2. Top 10 Suburbs By Average Price
#Create new dataframe - price and suburb columns
SuburbandPrice <- housing_data[c("Suburb","Price")]
#Calculate average price per suburb
top10avgprice <- SuburbandPrice %>% group_by(Suburb) %>%
  summarise(Average = sum(Price)/n()) %>%
  arrange(desc(Average)) %>% head(10)
#Plot
ggplot(top10avgprice, aes(reorder(Suburb, Average), Average, fill = Suburb))+
  geom_bar(stat = "identity")+
  theme(legend.position = "none")+
  labs(x = "Suburb", y = "Average Price Per Home",
       title = "Top 10 Suburbs by Average Price of Homes")

#3. Price Distribution of Homes
#Create new dataframe
DateandPrice <- housing_data[c("Suburb","Price","Date")] %>% na.omit()
#Plot
ggplot(DateandPrice, aes(Price))+
  geom_histogram(binwidth = 250000,color = "black", fill = "blue")+
  scale_x_continuous(breaks = c(1000000,2000000,3000000,4000000),
                     labels = c("$1m","$2m","$3m","$4m"))+
  ggtitle("Melbourne House Price Distribution")

#4. Distance from CBD Distribution
housing_data$Distance<- as.numeric(housing_data$Distance)
hist(housing_data$Distance, breaks = 40, xlim = c(0,50), ylim = c(0,3400),
     xlab = "Distance", col = "Blue", main = "Count of Distances from CBD", las =1)

#5. No. of Bedrooms Distribution
hist(housing_data$Bedroom2, breaks = 40, xlim = c(0,10), ylim = c(0,12000),
     xlab = "No. of Bedrooms", col = "grey", main = "Count of No. of Bedrooms", las =1)

#6. No. of Bathroom Distribution
hist(housing_data$Bathroom, breaks = 40, xlim = c(0,10), ylim = c(0,4000),
     xlab = "No. of Bathrooms", col = "Purple", main = "Count of No. of Bathrooms", las =1)

#7. No. of Car Lots Distribution
hist(housing_data$Car, breaks = 40, xlim = c(0,10), ylim = c(0,4000),
     xlab = "No. of Carslots", col = "Green", main = "Count of No. of Carslots", las =1)

#8. House Type Distribution
#Create New Dataframe
housedf2 <- housing_data %>%
  select(Price, Regionname,Type,Date)
#Plot
ggplot(housedf2, aes(Type, Price)) +
  geom_boxplot(outlier.colour = "blue") + 
  scale_x_discrete(labels = c('Houses','Townhouses','Units')) +
  scale_y_continuous(breaks=seq(0,10000000,1250000)) +
  xlab("Home Type") +
  ylab("Price") +
  ggtitle("Home Type Price Distribution")

#9. Price Distribution of Regions
#Create New Dataframe
housedf2 <- housing_data %>%
  select(Price, Regionname,Type,Date)
#Plot
ggplot(housedf2, aes(Regionname, Price)) +
  geom_boxplot(outlier.colour = "blue") +
  scale_y_continuous(breaks=seq(0,10000000,1250000)) +
  xlab("Region") +
  ylab("Price") +
  ggtitle("Price Distribution of Regions")

#10a. Correlation Matrix - All Homes
housedf3 <- housing_data %>%
  select(Rooms,Bedroom2,Bathroom,Distance,Car,Landsize,
         BuildingArea,YearBuilt,Propertycount,Price)
housedf3$Distance<- as.numeric(housedf3$Distance)
housedf3$Propertycount<- as.numeric(housedf3$Propertycount)
housedf3 <- na.omit(housedf3) 
corrplot(cor(as.matrix(housedf3)), method = "pie", type="lower")

#10b. Correlation Matrix - House Type Homes
housedfhome <- housing_data %>%filter(Type == "h") %>%
  select(Rooms,Bedroom2,Bathroom,Distance,Car,Landsize,
         BuildingArea,YearBuilt,Propertycount,Price)
housedfhome$Distance<- as.numeric(housedfhome$Distance)
housedfhome$Propertycount<- as.numeric(housedfhome$Propertycount)
housedfhome  <- na.omit(housedfhome ) 
corrplot(cor(as.matrix(housedfhome)), method = "pie", type="lower")

#10c. Correlation Matrix - Townhouse Type Homes
housedfhome <- housing_data %>%filter(Type == "t") %>%
  select(Rooms,Bedroom2,Bathroom,Distance,Car,Landsize,
         BuildingArea,YearBuilt,Propertycount,Price)
housedfhome$Distance<- as.numeric(housedfhome$Distance)
housedfhome$Propertycount<- as.numeric(housedfhome$Propertycount)
housedfhome  <- na.omit(housedfhome ) 
corrplot(cor(as.matrix(housedfhome)), method = "pie", type="lower")

#10d. Correlation Matrix - Unit Type Homes
housedfhome <- housing_data %>%filter(Type == "u") %>%
  select(Rooms,Bedroom2,Bathroom,Distance,Car,Landsize,
         BuildingArea,YearBuilt,Propertycount,Price)
housedfhome$Distance<- as.numeric(housedfhome$Distance)
housedfhome$Propertycount<- as.numeric(housedfhome$Propertycount)
housedfhome  <- na.omit(housedfhome ) 
corrplot(cor(as.matrix(housedfhome)), method = "pie", type="lower")
```

## Time Series

```{r, warning=FALSE, error=FALSE, message=FALSE}
datam<-housing_data

#change data type as.numeric(), as.character(), as.vector(), as.matrix(), as.data.frame)
#data pretreatment
dataT<-select(datam,Price,Date,Type)
dataT$Price<-as.numeric(dataT$Price)
dataT<-na.omit(dataT)

dataT1<-dataT %>% group_by(Date,Type) %>% summarise(Mean_sales=mean(Price))
#dataT1<-aggregate(dataT&Price,by=list(dataT&Date,dataT&Type),FUN=mean)
#conditional select mf[ mf$a == 16, ]

dataT2 <- subset(dataT1, Type == "u")
#View(dataT2)

dataT3<-select(dataT2,Date,Mean_sales)
#change date type
#dataT3$Date<-as.POSIXct(dataT3$Date,format = "%d/%m/%Y")
dataT3$Date<-as.Date(dataT3$Date, format = "%d/%m/%Y")
#View(dataT3)

dataT4<-xts(dataT3$Mean_sales, as.Date(dataT3$Date, format="%d/%m/%Y"))
#View(dataT4)

#If too many outliers, it could be changed to log format but in this case, it is not needed
#LAP<-log(dataT4)
AP<-dataT4
plot(dataT4)

#Decomposition of additive time series
#decomp<-decompose(AP)
#plot(decomp$figure,type='b',xlab='Month',ylab='Seasonality Index',col='blue',las=2)
#plot(decomp)
```

## Forecast

```{r, warning=FALSE, error=FALSE, message=FALSE}
#ARIMA- autoregressive integrated moving average
AP<-dataT4
model<-auto.arima(AP)
attributes(model)
model$coef

#ACF and PACF plots
#ACF plots display correlation between a series and its lag
acf(model$residuals,main='Correlogram')

#PACF plots display correlation between a series and its lag that explained by previous lag
pacf(model$residuals,main='Partial Correlogram')

#ljung-box test
Box.test(model$residuals, lag=20,type='Ljung-Box')

#residual plot
hist(model$residuals,col='red',xlab='Error',main='Histogram of Residuals',freq=FALSE)
lines(density(model$residuals))

#forecasr
f<-forecast(model,12)
autoplot(f)
accuracy(f)

#Forecasts using Exponential Smoothing  
skirtsseriesforecasts <- HoltWinters(AP, beta=FALSE,gamma=FALSE)
skirtsseriesforecasts$SSE
skirtsseriesforecasts2 <- forecast:::forecast.HoltWinters(skirtsseriesforecasts, h=12)
autoplot(skirtsseriesforecasts2)
accuracy(skirtsseriesforecasts2)

#Seasonal analysis 
#seasonm<-sarima(AP, 1,0,0,0,1,1,12)

#sarima.for(AP, 24, 0,1,1,0,1,1,12)

#Time series parameter
plot(AP, type="b")

#diff12 = diff(AP,12)
#acf2(diff12, 48)

#final
models = arima(AP, order = c(0,1,1), seasonal = list(order = c(0,1,1), period = 12))
predict(models, n.ahead=24) 
autoplot(forecast(models, 12))
accuracy(models)

#Extra code not related
#sarima(AP, 0,0,1,0,0,1,4)
#sarima(AP, 1,0,0,0,0,1,4)
#sarima(AP, 0,0,0,0,0,1,4)
#sarima(AP, 1,0,0,0,0,2,4)

#trunacate those under 0
themodel = arima(AP, order = c(1,0,0), seasonal = list(order = c(0,1,1), period = 12))
themodel
predict(themodel, n.ahead=12)
autoplot(forecast(themodel, n.ahead=12))

#monthly means to make the case that the data are seasonal.
flowm = matrix(AP, byrow=TRUE)
col.means=apply(flowm,2,mean)
plot(col.means,type="b", main="Monthly Means Plot for Flow", xlab="Month", ylab="Mean")

#seasonal analysis
models<-auto.arima(AP,seasonal=TRUE)
s<-forecast(models,12)
library(ggplot2)
autoplot(s)
accuracy(s)

#help(sarima)
models<-Arima(AP,order=c(0,1,2),seasonal=c(0,1,1))
s<-forecast(models,12)
library(ggplot2)
autoplot(s)
accuracy(s)

#fit<-stl(AP,s.wnidow="period")
#plot(fit)
#autoplot(s)
#accuracy(s)
```

## DL Modeling

#### Data Cleanup

```{r, warning=FALSE, error=FALSE, message=FALSE}
# Check for missing Price values and outliers
# Check for missing Price values
sum(is.na(housing_data$Price))

# Check for outliers
housing_data %>% filter(Price < 200000) %>% nrow()
housing_data %>% filter(Price < 300000) %>% nrow()
housing_data %>% filter(Price < 400000) %>% nrow()

housing_data %>% filter(Price > 3000000) %>% nrow()
housing_data %>% filter(Price > 4000000) %>% nrow()
housing_data %>% filter(Price > 5000000) %>% nrow()
housing_data %>% filter(Price > 6000000) %>% nrow()

# Since our target variable is Price, we drop observations without Price value.
# Start with full dataset
data_clean <- housing_data
dim(data_clean)

# Keep only records with non-null Price values
data_clean <- drop_na(data_clean, Price)
dim(data_clean)

# Keep only records with Price >= 300,000
data_clean <- data_clean %>% filter(Price >= 300000)
dim(data_clean)

# Keep only records with Price <= 5,000,000
data_clean <- data_clean %>% filter(Price <= 5000000)
dim(data_clean)
```

### Data Preparation

```{r, warning=FALSE, error=FALSE, message=FALSE}
# Divide the data into 10 buckets of equal sizes and assign Price Range value
num_buckets <- 10
data_clean$PriceRange <- as.numeric( cut_number( data_clean$Price, num_buckets) )

# Check the summary for the buckets
data_clean %>%
  group_by(PriceRange) %>%
  summarise(Count = n(), min = min(Price), max = max(Price), mean = mean(Price) )

str(data_clean)
```
### Features Selection

The following attributes were identified as having an impact on 
(are correlated to) the house price by our Data Analysis in section above.

(Type, Rooms, Bedroom2, Bathroom, Car, Landsize, YearBuilt, BuildingArea, Distance)

```{r, warning=FALSE, error=FALSE, message=FALSE}
# Data Cleanup for Independent attributes
data_numeric <- data_clean %>%
  select(Rooms, Bedroom2, Bathroom, Car, Landsize, YearBuilt, BuildingArea, Distance, Price, PriceRange)

# Convert non-numeric columns with numeric data to numeric.
data_numeric$Distance <- as.numeric(data_numeric$Distance)
str(data_numeric$Distance)

# Fill nulls value in numeric columns with column mean values
col_means <- lapply(data_numeric, mean, na.rm = TRUE)
datac_na  <- replace_na(data_numeric, col_means)

# Double-check null value in numeric columns
summarise_all(datac_na, funs( sum(is.na(.)) ) )

# Add categorical attribute (Type) to column selection
datac <- datac_na
datac$Type <- as.vector(data_clean$Type)

# Check distribution of categorical attribute (Type).
datac %>%
  group_by(Type) %>%
  summarise(Count = n(), min = min(Price), max = max(Price), mean = mean(Price) )

# Create dummies for categorical attribute (Type).
results <- fastDummies::dummy_cols(datac,select_columns='Type')

knitr::kable(results)

dataf <- results %>% select(Rooms, Bedroom2, Bathroom, Car, Landsize, YearBuilt, 
                            BuildingArea, Distance, Type_h,Type_t,Type_u, PriceRange)

dataf %<>% mutate_if(is.factor, as.numeric)

# Double-check null value in numeric columns
summarise_all(dataf, funs( sum(is.na(.)) ) )
```

### Model data

```{r, warning=FALSE, error=FALSE, message=FALSE}
model_data <- as.matrix(dataf)
dimnames(model_data) <- NULL

# Use 80/20 Train/Test Split
set.seed(123)
ind<-sample(2, nrow(model_data), replace=T, prob=c(0.8,0.2) )

train_set <- model_data[ind==1,1:8]
test_set  <- model_data[ind==2,1:8]

train_target <- model_data[ind==1,9]
test_target  <- model_data[ind==2,9]

#Normalization
mf <- colMeans(train_set)
sf <- apply(train_set, 2, sd)

train_set <- scale(train_set, center = mf, scale = sf)
test_set  <- scale(test_set, center = mf, scale = sf)

#Categorization
train_target = train_target -1
train_target_cat = to_categorical(train_target, num_classes = num_buckets)

test_target = test_target -1
test_target_cat  = to_categorical(test_target, num_classes = num_buckets)
```

### DL Model

Create a sequential model in Keras by stacking the layers sequentially.

- The model has 3 hidden layers.
- Keras builds an implicit input layer using the input_shape parameter.
- Last layer is the output layer.

```{r, warning=FALSE, error=FALSE, message=FALSE}

# Create a simple Neuralnet
nn <- neuralnet(PriceRange ~ Rooms+Bedroom2+Bathroom+Car+Landsize+YearBuilt+BuildingArea+Distance+Type_h+Type_t+Type_u,
               data = dataf,
               hidden = c(10,5),
               linear.output = F,
               lifesign = 'full',
               rep=1)
plot(nn,
     col.hidden = 'darkgreen',
     col.hidden.synapse = 'darkgreen',
     show.weights = F,
     information = F,
     fill = 'lightblue')

use_condaenv("r-tensorflow")

# Create a sequential model in Keras
model <- keras_model_sequential()
model %>% 
  layer_dense(units  = 100, activation = 'relu', input_shape = c(8) ) %>%
  layer_dense(units  = 50, activation = 'relu')  %>%
  layer_dense(units  = 20, activation = 'relu')  %>%
  layer_dense(units  = ncol(train_target_cat) )

# Take a look at the model and the shapes of the layers:
model

# Define the loss and optimizer functions and the metric to optimize.
compile(model, loss = 'categorical_crossentropy',
        optimizer = optimizer_rmsprop(), metrics = "accuracy")

# Fit (train) the model
model_hist <- fit(model, train_set, train_target_cat,
                  epochs=100, batch_size=32, validation_split=0.2)

# Plot fitting history
plot(model_hist)
```

### Evaluating the model

```{r, warning=FALSE, error=FALSE, message=FALSE}
# Using Testset data
evaluate(model, test_set, test_target_cat)

# Validate Model with unseen data
# Predict new cases
test_pred <- predict_classes(model, test_set)

pred_prob <- predict_proba(model, test_set)

# Confusion matrix
confusion_matrix = table(Predicted = test_pred, Actual = test_pred)
confusion_matrix

# Accuracy
accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)
accuracy
```