MULTIPLE LINEAR RGRESSION
Lumumba Wandera Victor
5/22/2022
#Data Importation
data<-read.csv("C:\\Users\\LUMUMBA\\Downloads\\training model.csv")
head(data,10)## year CPI Exch.Rate Lend.Int.Rates
## 1 1987 7.872727 16.45499 14.00000
## 2 1988 8.848083 17.74710 15.00000
## 3 1989 10.035029 20.57247 17.25000
## 4 1990 11.602322 22.91477 18.75000
## 5 1991 13.805882 27.50870 18.99750
## 6 1992 17.579356 32.21683 21.06750
## 7 1993 25.662147 58.00133 29.98917
## 8 1994 33.056538 56.05058 36.24000
## 9 1995 33.570345 51.42930 28.79583
## 10 1996 36.546050 57.11487 33.78667##attach(data)
attach(data)#View Dataset
head(data,10)## year CPI Exch.Rate Lend.Int.Rates
## 1 1987 7.872727 16.45499 14.00000
## 2 1988 8.848083 17.74710 15.00000
## 3 1989 10.035029 20.57247 17.25000
## 4 1990 11.602322 22.91477 18.75000
## 5 1991 13.805882 27.50870 18.99750
## 6 1992 17.579356 32.21683 21.06750
## 7 1993 25.662147 58.00133 29.98917
## 8 1994 33.056538 56.05058 36.24000
## 9 1995 33.570345 51.42930 28.79583
## 10 1996 36.546050 57.11487 33.78667tail(data,10)## year CPI Exch.Rate Lend.Int.Rates
## 17 2003 59.06105 75.93557 16.57333
## 18 2004 66.02725 79.17388 12.53167
## 19 2005 72.57203 75.55411 12.88250
## 20 2006 76.94992 72.10084 13.63533
## 21 2007 80.23609 67.31764 13.34034
## 22 2008 92.36287 69.17532 14.01694
## 23 2009 102.09550 77.35201 14.80454
## 24 2010 106.26492 79.23315 14.37150
## 25 2011 121.16535 88.81077 15.04676
## 26 2012 132.52856 84.52960 19.72341str(data)## 'data.frame': 26 obs. of 4 variables:
## $ year : int 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 ...
## $ CPI : num 7.87 8.85 10.04 11.6 13.81 ...
## $ Exch.Rate : num 16.5 17.7 20.6 22.9 27.5 ...
## $ Lend.Int.Rates: num 14 15 17.2 18.8 19 ...Plot the variables of interest
plot.ts(data$CPI)
plot.ts(data$Exch.Rate)
plot.ts(data$Lend.Int.Rates)
plot.ts(data$CPI,type="l")
#Activate important packages
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.8
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(vars)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: strucchange## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich##
## Attaching package: 'strucchange'## The following object is masked from 'package:stringr':
##
## boundary## Loading required package: urca## Loading required package: lmtestlibrary(olsrr)##
## Attaching package: 'olsrr'## The following object is masked from 'package:MASS':
##
## cement## The following object is masked from 'package:datasets':
##
## riverslibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somelibrary(ggplot2)
library(ggthemes)
library(grid)#Test stationarity
adf.test(data$CPI)## Warning in adf.test(data$CPI): p-value greater than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: data$CPI
## Dickey-Fuller = 1.5556, Lag order = 2, p-value = 0.99
## alternative hypothesis: stationaryadf.test(data$Exch.Rate)##
## Augmented Dickey-Fuller Test
##
## data: data$Exch.Rate
## Dickey-Fuller = -1.743, Lag order = 2, p-value = 0.6703
## alternative hypothesis: stationaryadf.test(data$Lend.Int.Rates)##
## Augmented Dickey-Fuller Test
##
## data: data$Lend.Int.Rates
## Dickey-Fuller = -2.0656, Lag order = 2, p-value = 0.5474
## alternative hypothesis: stationary#Difference the series
adf.test(diff(CPI))##
## Augmented Dickey-Fuller Test
##
## data: diff(CPI)
## Dickey-Fuller = -1.4862, Lag order = 2, p-value = 0.7681
## alternative hypothesis: stationaryCPI1<-diff(CPI)
head(CPI1,10)## [1] 0.9753552 1.1869462 1.5672929 2.2035602 3.7734740 8.0827912 7.3943910
## [8] 0.5138071 2.9757047 4.1376231Tst the ADF
adf.test(CPI1)##
## Augmented Dickey-Fuller Test
##
## data: CPI1
## Dickey-Fuller = -1.4862, Lag order = 2, p-value = 0.7681
## alternative hypothesis: stationaryadf.test(diff(CPI1))##
## Augmented Dickey-Fuller Test
##
## data: diff(CPI1)
## Dickey-Fuller = -3.1705, Lag order = 2, p-value = 0.1265
## alternative hypothesis: stationaryPlot and summarize the series
plot.ts(diff(CPI1))
adf.test(diff(CPI1))##
## Augmented Dickey-Fuller Test
##
## data: diff(CPI1)
## Dickey-Fuller = -3.1705, Lag order = 2, p-value = 0.1265
## alternative hypothesis: stationarySummarized the differenced series
summary(diff(CPI1))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -6.8806 -1.8506 0.2433 0.4328 1.7808 10.7310Exchange rate
adf.test(Exch.Rate)##
## Augmented Dickey-Fuller Test
##
## data: Exch.Rate
## Dickey-Fuller = -1.743, Lag order = 2, p-value = 0.6703
## alternative hypothesis: stationaryDifference the series
diff_Excha_Rate <-diff(Exch.Rate)Plot the series
plot.ts(diff_Excha_Rate)
#Data Visualization ### Consumer Price Index
hist(CPI)
hist(log(CPI))
hist(diff(CPI))
Exchange Rates
hist(Exch.Rate)
hist(log(Exch.Rate))
hist(diff(Exch.Rate))
Lending Interest Rates
hist(Lend.Int.Rates)
hist(log(Lend.Int.Rates))
hist(diff(Lend.Int.Rates))
#Test Normality
shapiro.test(CPI)##
## Shapiro-Wilk normality test
##
## data: CPI
## W = 0.93963, p-value = 0.1316shapiro.test(Exch.Rate)##
## Shapiro-Wilk normality test
##
## data: Exch.Rate
## W = 0.86288, p-value = 0.002555shapiro.test(Lend.Int.Rates)##
## Shapiro-Wilk normality test
##
## data: Lend.Int.Rates
## W = 0.86773, p-value = 0.003203#ESTIMATE REGRESSION EQUATION
model<-lm(log(CPI)~log(Exch.Rate)+log(Lend.Int.Rates),data=data)
model##
## Call:
## lm(formula = log(CPI) ~ log(Exch.Rate) + log(Lend.Int.Rates),
## data = data)
##
## Coefficients:
## (Intercept) log(Exch.Rate) log(Lend.Int.Rates)
## -0.7043 1.5149 -0.5520Model Summary
summary(model)##
## Call:
## lm(formula = log(CPI) ~ log(Exch.Rate) + log(Lend.Int.Rates),
## data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3328 -0.1943 0.0154 0.1211 0.5152
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.70426 0.52767 -1.335 0.195053
## log(Exch.Rate) 1.51494 0.08531 17.758 6.35e-15 ***
## log(Lend.Int.Rates) -0.55204 0.13754 -4.014 0.000544 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.223 on 23 degrees of freedom
## Multiple R-squared: 0.9348, Adjusted R-squared: 0.9291
## F-statistic: 164.8 on 2 and 23 DF, p-value: 2.317e-14Activate the library stargazer
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerstargazer(model,type="text")##
## ===============================================
## Dependent variable:
## ---------------------------
## log(CPI)
## -----------------------------------------------
## log(Exch.Rate) 1.515***
## (0.085)
##
## log(Lend.Int.Rates) -0.552***
## (0.138)
##
## Constant -0.704
## (0.528)
##
## -----------------------------------------------
## Observations 26
## R2 0.935
## Adjusted R2 0.929
## Residual Std. Error 0.223 (df = 23)
## F Statistic 164.844*** (df = 2; 23)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01#Model Diagnostic
1. Test for the presence of outliers in the model
outlierTest(model)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 26 2.702426 0.013008 0.3382The p-value for Boniferron test statistics shows that there are outliers in the data set
qqPlot(model, main= "QQ Plot Showing the Possible Presence of outliers")
## [1] 18 26leveragePlots(model)
2. test the presence of multicollinearity in the model
vif(model)## log(Exch.Rate) log(Lend.Int.Rates)
## 1.000162 1.000162VIF of 1.000162 is an indication that predictors are no correlated
Testing for the presence of autocorrelation
3. Testing for heteroscedasticty
ncvTest(model)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 5.509929, Df = 1, p = 0.018909the results show that the variance of the error terms is not constant
Spread Level Plot
spreadLevelPlot(model)
##
## Suggested power transformation: -2.486722Dwattson test statistics
durbinWatsonTest(model)## lag Autocorrelation D-W Statistic p-value
## 1 0.6176191 0.5324069 0
## Alternative hypothesis: rho != 0The results shows that there is a correlation of the regression residuals
Correcting Econometric problems
#Heteroscedasticity
coeftest(model, hccm(model, type = "hc0"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.704260 0.315018 -2.2356 0.03537 *
## log(Exch.Rate) 1.514945 0.057345 26.4181 < 2.2e-16 ***
## log(Lend.Int.Rates) -0.552044 0.109898 -5.0232 4.403e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Present the results usng stargazer
stargazer(coeftest(model, hccm(model, type = "hc0")),type="text")##
## ===============================================
## Dependent variable:
## ---------------------------
##
## -----------------------------------------------
## log(Exch.Rate) 1.515***
## (0.057)
##
## log(Lend.Int.Rates) -0.552***
## (0.110)
##
## Constant -0.704**
## (0.315)
##
## ===============================================
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01VISUAL REPRESENTATION OF THE TIME SERIES DATA USING GGPLOT2 AND GGTHEMES
#We need to declare our data be a time series (Consumer Price Index)
CPI<-ts(data$CPI,start=1987,frequency = 1)Plot the series
plot.ts(CPI,type="l",main="Time Series plot CPI",xlab="Year",ylab="Consumer Price Index")
Plot the differenced series
plot.ts(diff(CPI),type="l",main="Time Series plot CPI",xlab="Year",ylab="Consumer Price Index")
Exchange Rate
Declare the series time series
Exch.Rate<-ts(data$Exch.Rate,start=1987,frequency = 1)Plot the series
plot.ts(Exch.Rate,type="l",main="Time Series plot Exch.Rate",xlab="Year",ylab="Exchange Rate")
plot.ts(diff(Exch.Rate),type="l",main="Time Series plot Exch.Rate",xlab="Year",ylab="Exchange Rate")
Lending Interest Rates
Lend.Int.Rates<-ts(data$Lend.Int.Rates,start=1987,frequency = 1)Plot the series
plot.ts(Lend.Int.Rates,type="l",main="Time Series plot Lend.Int.Rates",xlab="Year",ylab="Lend.Int.Rates")
plot.ts(diff(Lend.Int.Rates),type="l",main="Time Series plot Lend.Int.Rates",xlab="Year",ylab="Lend.Int.Rates")
DATA VISUALIZATION USING GGPLOT2 AND GGTHEMES
#####1
ggplot(data=data,aes(x=year,y=CPI))+geom_line()
#####2
ggplot(data=data,aes(x=year,y=CPI))+geom_line()+
labs(title="Time Series plot of CPI",
caption="source:World Bank",
y="Consumer Price Index", x="Year",
color=3) + # title and caption
theme(axis.text.x = element_text(angle = 0, vjust=0.5, size = 12), # rotate x axis text
axis.title=element_text(size=12,face="bold"),
panel.grid = element_blank())+
#theme(panel.grid.minor = element_blank())+#turn off minor grid(to run remove #be4 theme)
theme(legend.text = element_text(size=12,face="bold"))+
theme_set(theme_economist())
#####3
ggplot(data=data,aes(x=year,y=CPI))+geom_line()+
labs(title="Time Series plot of CPI",
caption="source:World Bank 2018", y="Consumer Price Index", x="Year")
4
ggplot(data=data,aes(x=year,y=Exch.Rate))+geom_line()+
labs(title="Time Series plot of Exhange Rate",
caption="source:World Bank 2018", y="Exhange Rate", x="Year")
5
ggplot(data=data,aes(x=year,y=Lend.Int.Rates))+geom_line()+
labs(title="Time Series plot of Lend.Int.Rates",
caption="source:World Bank 2018", y="Lend.Int.Rates", x="Year")
Rename time variable
date<-seq(as.Date("1987-01-01"),by="1 year",length.out=length(data$year))
date## [1] "1987-01-01" "1988-01-01" "1989-01-01" "1990-01-01" "1991-01-01"
## [6] "1992-01-01" "1993-01-01" "1994-01-01" "1995-01-01" "1996-01-01"
## [11] "1997-01-01" "1998-01-01" "1999-01-01" "2000-01-01" "2001-01-01"
## [16] "2002-01-01" "2003-01-01" "2004-01-01" "2005-01-01" "2006-01-01"
## [21] "2007-01-01" "2008-01-01" "2009-01-01" "2010-01-01" "2011-01-01"
## [26] "2012-01-01"Plot the series
ggplot(data=data,aes(x=date))+
geom_line(aes(y=Exch.Rate,colour="Exhange Rate"))+
geom_line(aes(y=Lend.Int.Rates,colour="Lending Interest Rates"))+
geom_line(aes(y=CPI,colour="Consumer Price Index"))+
labs(title="Trends of CPI,Interest Rates and Exchange Rates",
caption="", y="Rate", x="Time in Years", color="Key")+
scale_x_date( date_labels = "%Y", breaks = "1 year")+
theme(axis.text.x = element_text(angle = 90, vjust=0.5, size = 8))
theme_set(theme_dark())
theme_set(theme_economist())
theme_set(theme_base())Create plot1
p1<-ggplot(data=data,aes(x=date,y=CPI))+geom_line()+
labs(title="Consumer Price Index",
caption="", y="Consumer Price Index", x="Time in Years", color=3)+
scale_x_date( date_labels = "%Y-%b", breaks = "1 years")+
theme(axis.text.x = element_text(angle = 90, vjust=0.5, size = 8))Create plot2
p2<-ggplot(data=data,aes(x=date,y=Exch.Rate))+geom_line()+
labs(title="Exchange Rate",
caption="", y="Exchange Rate", x="Time in Years", color="Key")+
scale_x_date( date_labels = "%Y-%b", breaks = "1 years")+
theme(axis.text.x = element_text(angle = 90, vjust=0.5, size = 8))Create Plot3
p3<-ggplot(data=data,aes(x=date,y=Lend.Int.Rates))+geom_line()+
labs(title="Lending Interest Rates",
caption="", y="Lending Interest Rates", x="Time in Years", color="Key")+
scale_x_date( date_labels = "%Y-%b", breaks = "1 years")+
theme(axis.text.x = element_text(angle = 90, vjust=0.5, size = 8))Plot the three plots
grid.newpage()
grid.draw(rbind(ggplotGrob(p1),ggplotGrob(p2),ggplotGrob(p3),size="last"))
Plot the 1st and 3rd plots
grid.newpage()
grid.draw(rbind(ggplotGrob(p1),ggplotGrob(p3),size="last"))
Plot the 2nd and 3rd plots
grid.newpage()
grid.draw(rbind(ggplotGrob(p2),ggplotGrob(p3),size="last"))
Plot the 1st plot
grid.newpage()
grid.draw(rbind(ggplotGrob(p1),size="last"))
Plot the 2nd plot
grid.newpage()
grid.draw(rbind(ggplotGrob(p2),size="last"))
Plot the third Plot
grid.newpage()
grid.draw(rbind(ggplotGrob(p3),size="last"))