Car Sales and Macro-Economics
Christophe Duong
12/21/2014
1 Data Source and Description
Two datasets were extracted manually from knoema.com:
World Motor Vehicle Sales by country and type 2005-2012 (Publication date: 25 February, 2013): http://knoema.com/lnubnsd/world-motor-vehicle-sales-by-country-and-type-2005-2012
The World Economic Outlook (WEO) database contains selected macroeconomic data series from the statistical appendix of the World Economic Outlook report, which presents the IMF staff’s analysis and projections of economic developments at the global level, in major country groups and in many individual countries: (Publication date: 07 October, 2014) http://knoema.com/IMFWEO2014Oct/imf-world-economic-outlook-october-2014
After retrieving them, these datasets were merged and combined into a single csv file that we are now using in the following analysis using R where we try to build a Linear Regression model using OLS to predict Car Sales figures depending on macroeconomics measures of any country in the world.
2 Basic Analysis
2.1 Data loading and preparation
Load the data, shorten the column names to get cleaner output and remove outliers from all country (named under ‘world’):
df <- read.csv('data.csv')
names(df)[names(df) == 'Gross.domestic.product..current.prices..U.S..dollars.'] <- 'GDP.USD'
names(df)[names(df) == 'Current.account.balance..Percent.of.GDP.'] <- 'Current.Acc.Bal.PercGDP'
names(df)[names(df) == 'Current.account.balance..U.S..dollars.'] <- 'Current.Acc.Bal.USD'
names(df)[names(df) == 'Employment..Persons.'] <- 'Employment'
names(df)[names(df) == 'Export.volume.of.goods.and.services..Percent.change.'] <- 'Export'
names(df)[names(df) == 'Gross.domestic.product.based.on.purchasing.power.parity..PPP..valuation.of.country.GDP..Current.international.dollar.'] <- 'GDP.PPP'
names(df)[names(df) == 'Gross.domestic.product.per.capita..current.prices..U.S..dollars.'] <- 'GDP.per.capita'
names(df)[names(df) == 'Gross.domestic.product..constant.prices..Percent.change.'] <- 'GDP.PercChange'
names(df)[names(df) == 'Gross.domestic.product..deflator..Index.'] <- 'GDP.deflatorIndex'
names(df)[names(df) == 'Gross.national.savings..Percent.of.GDP.'] <- 'National.Savings'
names(df)[names(df) == 'Import.volume.of.goods.and.services..Percent.change.'] <- 'Import'
names(df)[names(df) == 'Inflation..average.consumer.prices..Index.'] <- 'Inflation.Index'
names(df)[names(df) == 'Inflation..average.consumer.prices..Percent.change.'] <- 'Inflation.PercChange'
names(df)[names(df) == 'Investment..Percent.of.GDP.'] <- 'Investment'
names(df)[names(df) == 'Output.gap.in.percent.of.potential.GDP..Percent.of.potential.GDP.'] <- 'Output.Gap'
names(df)[names(df) == 'Population..Persons.'] <- 'Population'
names(df)[names(df) == 'Trade.volume.of.goods.and.services..Percent.change.'] <- 'Trade.Volume'
names(df)[names(df) == 'Unemployment.rate..Percent.of.total.labor.force.'] <- 'Unemployment.Rate'
df <- df[df$GDP.USD < 40000,]
todrop <- c("Commercial.vehicles","Passengers.Cars", "Trade.Volume")
df <- df[,!(names(df) %in% todrop)]2.2 Basic Scatter Plot
We want to check if some macroeconomics metrics would be linearly dependent with car sales in the world. Let’s first look at a scatter plot:
library(ggplot2)
qplot(df$All.vehicles, df$GDP.USD) +
xlab("Amount of All Vehicles Sold") +
ylab("GDP in US Dollars")
2.3 Fit a Simple Linear Regression Model
Let’s fit a linear regression model through that data and look at the resulting residuals.
As we are using only one independent variable, it is called a simple linear regression:
excludedCountries <- c('japan', 'china', 'united states of america');
dfExcluded <- df[!(df$country %in% excludedCountries),]
dfNonExcluded <- df[(df$country %in% excludedCountries),]
fit <- lm(data=dfExcluded, All.vehicles ~ GDP.USD)
summary(fit)##
## Call:
## lm(formula = All.vehicles ~ GDP.USD, data = dfExcluded)
##
## Residuals:
## Min 1Q Median 3Q Max
## -739231 -43604 -778 10171 1479081
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9640.85 7644.06 -1.261 0.208
## GDP.USD 1133.73 12.62 89.804 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 220400 on 1030 degrees of freedom
## Multiple R-squared: 0.8867, Adjusted R-squared: 0.8866
## F-statistic: 8065 on 1 and 1030 DF, p-value: < 2.2e-16yhat <- predict(fit)
yhat2 <- predict(fit, newdata = dfNonExcluded)We actually identified visually a certain number of countries such as Japan, China and USA to behave quite differently from the other countries of the world as their car sales are not as linearly dependent as the others. Therefore, we choose to display the residual bars in green for all the countries other than these ‘outlier’ countries which are displayed in red:
qplot(x = GDP.USD, y = All.vehicles, data=df) +
geom_smooth(method="lm", se=FALSE, size=1, data=df) +
geom_segment(data= dfExcluded, x=dfExcluded$GDP.USD,
y=dfExcluded$All.vehicles,
xend=dfExcluded$GDP.USD,
yend=yhat, colour=I("green"),alpha=0.5) +
geom_point(data=dfNonExcluded, shape=23) +
geom_segment(data= dfNonExcluded, x=dfNonExcluded$GDP.USD,
y=dfNonExcluded$All.vehicles,
xend=dfNonExcluded$GDP.USD,
yend=yhat2, colour=I("red"),alpha=0.5) +
ylab("Amount of All Vehicles Sold") +
xlab("GDP in US Dollars")
Using some various built-in regression diagnostics plots:
par(mfrow=c(2,2))
plot(fit)
We can already notice that some of the assumptions on this model based on linear regression are not respected (normality for example)
3 Data Analysis using more variables
3.1 Fit a Multiple Linear Regression Model
Now, let’s consider all the independent variables together, therefore we are running a multiple linear regression instead:
fitMore <- lm(data=df, All.vehicles ~ . - year - X - country)
summary(fitMore)##
## Call:
## lm(formula = All.vehicles ~ . - year - X - country, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2231848 -165246 -14508 128200 1861277
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.755e+05 4.828e+05 1.606 0.10988
## Current.Acc.Bal.PercGDP -3.512e+05 1.194e+05 -2.942 0.00367 **
## Current.Acc.Bal.USD -3.959e+03 4.603e+02 -8.601 2.92e-15 ***
## Employment 1.449e+05 1.900e+04 7.624 1.14e-12 ***
## Export -1.030e+04 8.223e+03 -1.253 0.21188
## GDP.PPP 2.463e+02 2.013e+02 1.223 0.22275
## GDP.per.capita -4.718e+00 2.268e+00 -2.080 0.03886 *
## GDP.PercChange 1.117e+03 1.682e+04 0.066 0.94712
## GDP.USD -3.894e+02 1.895e+02 -2.055 0.04124 *
## GDP.deflatorIndex 4.590e+03 3.111e+03 1.475 0.14183
## National.Savings 3.818e+05 1.193e+05 3.199 0.00161 **
## Import 1.271e+04 7.820e+03 1.625 0.10576
## Inflation.Index -7.296e+03 4.673e+03 -1.561 0.12014
## Inflation.PercChange 1.967e+04 2.641e+04 0.745 0.45727
## Investment -3.942e+05 1.191e+05 -3.311 0.00111 **
## Output.Gap 2.422e+04 1.317e+04 1.840 0.06740 .
## Population -1.790e+04 1.089e+04 -1.643 0.10200
## Unemployment.Rate -8.717e+03 1.197e+04 -0.728 0.46748
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 402000 on 190 degrees of freedom
## (848 observations deleted due to missingness)
## Multiple R-squared: 0.9826, Adjusted R-squared: 0.9811
## F-statistic: 632.5 on 17 and 190 DF, p-value: < 2.2e-163.2 Stepwise Linear Regression Model
The previous model was ending up using a lot of independent variables that are not significant. Let’s optimize the model by trying to add or remove some of them and optimize by looking at the Akaike information criterion (AIC) value which tries to deal with the concepts of goodness of fit of the model and its complexity:
library(MASS)
fitStep <- stepAIC(fitMore, direction="both", trace=0)
summary(fitStep)##
## Call:
## lm(formula = All.vehicles ~ Current.Acc.Bal.PercGDP + Current.Acc.Bal.USD +
## Employment + GDP.per.capita + GDP.USD + National.Savings +
## Import + Investment + Output.Gap + Unemployment.Rate, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2336174 -150819 -15313 142328 1862236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.283e+05 2.445e+05 2.161 0.031919 *
## Current.Acc.Bal.PercGDP -3.818e+05 1.156e+05 -3.302 0.001139 **
## Current.Acc.Bal.USD -3.937e+03 3.659e+02 -10.760 < 2e-16 ***
## Employment 1.185e+05 7.956e+03 14.891 < 2e-16 ***
## GDP.per.capita -4.534e+00 2.154e+00 -2.105 0.036587 *
## GDP.USD -3.165e+02 8.454e+01 -3.744 0.000238 ***
## National.Savings 4.116e+05 1.156e+05 3.562 0.000461 ***
## Import 5.422e+03 3.401e+03 1.594 0.112520
## Investment -4.206e+05 1.154e+05 -3.646 0.000341 ***
## Output.Gap 2.292e+04 1.146e+04 2.000 0.046885 *
## Unemployment.Rate -2.345e+04 9.434e+03 -2.486 0.013761 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 400300 on 197 degrees of freedom
## (848 observations deleted due to missingness)
## Multiple R-squared: 0.9821, Adjusted R-squared: 0.9812
## F-statistic: 1084 on 10 and 197 DF, p-value: < 2.2e-16We can now notice that the R-Squared and model’s overal p-value has been barely affected but we went from a F-Statistic of 632.5 to 1084 with the stepwise approach.
As a test, let’s try to remove even some more of the non-significant variables and see the effect on the R-Squared value:
fitStep <- lm(data=df, All.vehicles ~ Current.Acc.Bal.PercGDP + Current.Acc.Bal.USD + Employment + GDP.USD + National.Savings + Investment + Output.Gap)
summary(fitStep)##
## Call:
## lm(formula = All.vehicles ~ Current.Acc.Bal.PercGDP + Current.Acc.Bal.USD +
## Employment + GDP.USD + National.Savings + Investment + Output.Gap,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2447926 -134192 -10383 140473 1921686
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39031.6 194165.8 0.201 0.840886
## Current.Acc.Bal.PercGDP -364888.7 116672.4 -3.127 0.002026 **
## Current.Acc.Bal.USD -3703.1 362.8 -10.207 < 2e-16 ***
## Employment 120797.7 7913.9 15.264 < 2e-16 ***
## GDP.USD -329.4 84.4 -3.903 0.000130 ***
## National.Savings 390335.5 116310.6 3.356 0.000946 ***
## Investment -393690.3 115914.1 -3.396 0.000824 ***
## Output.Gap 34058.5 10625.2 3.205 0.001570 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 408300 on 200 degrees of freedom
## (848 observations deleted due to missingness)
## Multiple R-squared: 0.9811, Adjusted R-squared: 0.9805
## F-statistic: 1487 on 7 and 200 DF, p-value: < 2.2e-16The Adjusted R-Squared is now slightly decreasing from 0.9812 to 0.9805 but the F-Statistic is again increasing from 1084 to 1487.
3.3 Eliminating correlation and looking at interactions
By studying the VIF between the variables, we can pursue further and decide to remove some of the correlated variables too and then study interactions with the remaining ones:
fitInteraction <- lm(data=df, All.vehicles ~ GDP.USD*Investment + Output.Gap)
summary(fitInteraction)##
## Call:
## lm(formula = All.vehicles ~ GDP.USD * Investment + Output.Gap,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1416561 -199882 -30574 151027 1952661
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.155e+06 1.978e+05 5.839 2.06e-08 ***
## GDP.USD -4.345e+02 8.372e+01 -5.190 5.08e-07 ***
## Investment -5.554e+04 8.958e+03 -6.200 3.10e-09 ***
## Output.Gap 4.841e+04 1.026e+04 4.719 4.40e-06 ***
## GDP.USD:Investment 6.997e+01 4.089e+00 17.113 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 403900 on 203 degrees of freedom
## (848 observations deleted due to missingness)
## Multiple R-squared: 0.9813, Adjusted R-squared: 0.9809
## F-statistic: 2659 on 4 and 203 DF, p-value: < 2.2e-16library(effects)
plot(effect("GDP.USD:Investment", fitInteraction, xlevels=list(wt=c(2,3))), multiline=TRUE)
fitStep <- fitInteractionAgain the R-Squared is not significantly affected but the F-Statistic gain a bigger amount of increase up to 2659. We can also conclude from the interaction graph that depending on the GDP level, when Investments made by the goverment are higher, there is also higher car sales in that country. This effect is more pronounced when the GDP is higher than when it is lower.
3.4 Regression Diagnostics
3.4.1 Checking Linear Relation
If the relation between the IVs and the DV is linear, no pattern should appear in the scatter plot:
library(car)##
## Attaching package: 'car'
##
## The following object is masked from 'package:effects':
##
## Prestigeqplot(predict(fitStep), resid(fitStep), geom="point") + geom_hline(yintercept=0)
As it seems to exist a slight pattern, the linear relation in this data set may not hold true.
3.4.2 Checking Normality
x <- resid(fitStep)
qqPlot(fitStep,main = "QQ-Plot")
qplot(x,geom="blank") +
geom_histogram( colour=I("white"), aes(y=..density..)) +
stat_function(fun=dnorm, aes(colour="Normal"),arg=list(mean=mean(resid(fitStep)),
sd=sd(resid(fitStep))))## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
It almost look normal but because of some outliers as denoted previously in the scatter plot, normality may not be true either for all the countries.
3.4.3 Checking Homoscedasticity
library(lmtest)## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericqplot(predict(fitStep),resid(fitStep), geom="point")
spreadLevelPlot(fitStep)## Warning in spreadLevelPlot.lm(fitStep): 24 negative fitted values removed
##
## Suggested power transformation: 0.4082004bptest(fitStep)##
## studentized Breusch-Pagan test
##
## data: fitStep
## BP = 48.7869, df = 4, p-value = 6.468e-10The line doesn’t seem to be horizontal and the statistic test of Breusch–Pagan seems to indicate heteroskedasticity in our model with a significant p-value below our alpha threshold so we may not satisfy the homogeneity of variance assumption either.
3.4.4 Checking Independence
durbinWatsonTest(fitStep)## lag Autocorrelation D-W Statistic p-value
## 1 0.4511046 1.070736 0
## Alternative hypothesis: rho != 0The Durbin-Watson test value of 1.367461 is not close to 2 but the p-value is below alpha so we can conclude that errors are not correlated. (as the value is not less than 1.0, there may not be cause for alarm)
4 Conclusion
As a conclusion, even though we may have obtained an amazing Ordinary Least Square Linear Regression model with an adjusted R-Squared of 0.9809, we failed to verify multiple of its assumptions so the built model may not be usable to predict and plan marketing strategies of Car Sales per country based on macroeconomics measures.