HUL715 Paper Replication
Ayanik Anwesh Patra
2022-19-02
Does Stock Market Development Cause Economic Growth?
- For the purpose of the study we take inspiration from the paper by Dev and Mukherjee
library(readxl)
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(aod)
library(pdfetch)
library(lubridate)## Loading required package: timechange##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(vars)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: strucchange## Loading required package: sandwich## Loading required package: urca## Loading required package: lmtestlibrary(FinAna)##
## Attaching package: 'FinAna'## The following object is masked from 'package:dplyr':
##
## desclibrary(strucchange)
library(sjPlot)directory <- "/cloud/project/data"DATA CONCERN
- Data on GDP is taken from Fred St. Louis database
- Data on Market capitalization and turnover (trading volume) are taken from CMIE PROwess IQ
- Data on Stock Market Volatility (S&P BSE Sensex) is taken from Yahoo_finance & Google Finance
gdpind <- as.data.frame(read.csv(file.path(directory,"NAEXKP01INQ652S.csv"))) # Real GDP
mktcap <- as.data.frame(read_excel(file.path(directory,"Final.xlsx"),sheet = "Market Capitalisaton"))
trnvr <- as.data.frame(read_excel(file.path(directory,"Final.xlsx"),sheet = "Turnover"))# Real GDP
gdpind$qy=as.yearqtr(as.Date( gdpind$DATE, "%Y-%m-%d"))
gdpind.zoo <- zoo(gdpind$NAEXKP01INQ652S,order.by = gdpind$qy)
#Market Capitalization
mktcap$qy=as.yearqtr(as.Date( mktcap$DATE, "%Y-%m-%d"))
mktcap.zoo <- zoo(mktcap$MKTCAP,order.by = mktcap$qy)
#Turnover/ Trading Volume
trnvr$qy=as.yearqtr(as.Date( trnvr$DATE, "%Y-%m-%d"))
trnvr.zoo <- zoo(trnvr$TRNR,order.by = trnvr$qy)#Stock Volatility
tickers <- c("^BSESN")
volat <- pdfetch_YAHOO(tickers,fields = "close",from="1996-01-01", to="2021-12-31")
write.zoo(x=volat,file = "Volatillity.csv",col.names = c("DATE","BESN"))
volatility <- as.data.frame(read.csv(file.path(directory,"Volatility.csv")))
data <- volatility
data$year <- strftime(data$DATE,"%Y")
data$month <- strftime(data$DATE,"%m")
data_aggr <- aggregate(BESN~month+year,data,FUN = mean) # Finding the mean of daily data for using as monthly average
data_aggr$Date <- as.yearmon(paste(data_aggr$year, data_aggr$month), "%Y %m ")
data_aggr$qy=as.yearqtr(as.yearmon(data_aggr$Date,"%m %Y"))
data_aggr.zoo <- zoo(data_aggr$BESN,order.by = data_aggr$qy)## Warning in zoo(data_aggr$BESN, order.by = data_aggr$qy): some methods for "zoo"
## objects do not work if the index entries in 'order.by' are not unique# Stock Market Volatility(SMV)
smv <- data_aggr
smv <- aggregate(BESN~qy,data_aggr,FUN = mean)
smv.zoo <- zoo(smv$BESN,order.by = smv$qy)par(mfrow=c(2,2))
plot(gdpind.zoo,main="Real GDP",ylab="Real GDP",xlab="Year")
plot(mktcap.zoo,main="Market Capitalisation",ylab="Mkt Cap",xlab="Year")
plot(trnvr.zoo,main="Market Turnover",ylab="Turnover",xlab="Year")
plot(smv.zoo,main="Stock market Volatility",ylab="Volatility",xlab="Year")
par(mfrow=c(2,2))
acf(ts(gdpind.zoo, f=1),lag.max = 32, main="ACF of GDP",
xaxt="n", xlab="Lag (quarter)", ylim=c(-.2,1))
axis(1, at=seq(0, 32, 1))
acf(ts(mktcap.zoo, f=1),lag.max = 32, main="ACF of Market Capitalisation",
xaxt="n", xlab="Lag (quarter)", ylim=c(-.2,1))
axis(1, at=seq(0, 32, 1))
acf(ts(trnvr.zoo, f=1),lag.max = 32, main="ACF of Turnover",
xaxt="n", xlab="Lag (quarter)", ylim=c(-.2,1))
axis(1, at=seq(0, 32, 1))
acf(ts(smv.zoo, f=1),lag.max = 32, main="ACF of Volatility",
xaxt="n", xlab="Lag (quarter)", ylim=c(-.2,1))
axis(1, at=seq(0, 32, 1))
mtext("My Multiplot Title", # Add main title
side = 5,
line = -2,
outer = TRUE)
- Clearly this indicates all the variables are stochastic in nature
Checking for structural breaks
str.gdp <- ts(gdpind$NAEXKP01INQ652S,start=1996,frequency = 4)
dfgdp <- breakpoints(str.gdp~1)
summary(dfgdp) # potential breakpoints of real GDP##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = str.gdp ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 62
## m = 2 44 76
## m = 3 36 57 79
## m = 4 35 54 70 85
## m = 5 26 41 56 71 86
##
## Corresponding to breakdates:
##
## m = 1 2011(2)
## m = 2 2006(4) 2014(4)
## m = 3 2004(4) 2010(1) 2015(3)
## m = 4 2004(3) 2009(2) 2013(2) 2017(1)
## m = 5 2002(2) 2006(1) 2009(4) 2013(3) 2017(2)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 9.956e+27 2.237e+27 8.718e+26 4.930e+26 3.548e+26 2.864e+26
## BIC 6.588e+03 6.440e+03 6.351e+03 6.300e+03 6.275e+03 6.262e+03str.mcap <- ts(mktcap$MKTCAP,start = 1996,frequency = 4)
dfmcap <- breakpoints(str.mcap~1)
summary(dfmcap) # potential breakpoints of market capitalization##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = str.mcap ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 73
## m = 2 45 84
## m = 3 43 72 87
## m = 4 38 54 72 87
## m = 5 24 39 54 72 87
##
## Corresponding to breakdates:
##
## m = 1 2014(1)
## m = 2 2007(1) 2016(4)
## m = 3 2006(3) 2013(4) 2017(3)
## m = 4 2005(2) 2009(2) 2013(4) 2017(3)
## m = 5 2001(4) 2005(3) 2009(2) 2013(4) 2017(3)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 3.824e+15 1.121e+15 5.748e+14 4.131e+14 3.779e+14 3.738e+14
## BIC 3.553e+03 3.435e+03 3.374e+03 3.349e+03 3.349e+03 3.358e+03str.trnvr <- ts(trnvr$TRNR,start = 1996,frequency = 4)
dftrnvr <- breakpoints(str.trnvr~1)
summary(dftrnvr) # potential breakpoints of turnover or total traded value##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = str.trnvr ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 86
## m = 2 45 87
## m = 3 40 72 87
## m = 4 30 45 72 87
## m = 5 15 30 45 72 87
##
## Corresponding to breakdates:
##
## m = 1 2017(2)
## m = 2 2007(1) 2017(3)
## m = 3 2005(4) 2013(4) 2017(3)
## m = 4 2003(2) 2007(1) 2013(4) 2017(3)
## m = 5 1999(3) 2003(2) 2007(1) 2013(4) 2017(3)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 2.582e+10 8.991e+09 4.798e+09 4.472e+09 4.238e+09 4.192e+09
## BIC 2.315e+03 2.214e+03 2.158e+03 2.160e+03 2.164e+03 2.172e+03str.smv <- ts(smv$BESN,start = 1997,frequency = 4)
dfsmv <- breakpoints(str.smv~1)
summary(dfsmv) # potential breakpoints of stock market volatility##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = str.smv ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 67
## m = 2 37 79
## m = 3 35 67 84
## m = 4 34 49 67 84
## m = 5 20 34 49 67 84
##
## Corresponding to breakdates:
##
## m = 1 2013(3)
## m = 2 2006(1) 2016(3)
## m = 3 2005(3) 2013(3) 2017(4)
## m = 4 2005(2) 2009(1) 2013(3) 2017(4)
## m = 5 2001(4) 2005(2) 2009(1) 2013(3) 2017(4)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 1.797e+10 5.221e+09 2.516e+09 1.442e+09 1.254e+09 1.238e+09
## BIC 2.152e+03 2.040e+03 1.978e+03 1.932e+03 1.928e+03 1.936e+03Plotting Structural Breaks (possible breakpoints)
par(mfrow=c(2,2))
ci.dfgdp <- confint(dfgdp)
ts.plot(str.gdp) #GDP
lines(ci.dfgdp)
ci.dfmcap <- confint(dfmcap)
ts.plot(str.mcap) #Market Capitalisation
lines(ci.dfmcap)
ci.dftrnvr <- confint(dftrnvr)
ts.plot(str.trnvr) #Turnover
lines(ci.dftrnvr)
ci.dfsmv <- confint(dfsmv)
ts.plot(str.smv) #Stock Volatility
lines(ci.dfsmv)
mtext("STRUCTURAL BREAKS PLOTS", side = 3, line = -1, outer = TRUE)
* Clearly there is presence of structural breaks in the variables , so
ordinary ADF (Augmented Dickey Fuller) tests are biased.
#Checking Unit root
- How to check Stationarity
- Ho:The data has unit root or the data is not stationary.
- Perron showed that in presence of structural breaks ADF test doesn’t distinguish between stochastic and deterministic trends. -When structural breaks are present ADF tests are biased to non-rejection of unit root
# We'll use will use Phillips Perron Test as it is not susceptible to structural breaks like ADF test
PP.test(gdpind$NAEXKP01INQ652S)##
## Phillips-Perron Unit Root Test
##
## data: gdpind$NAEXKP01INQ652S
## Dickey-Fuller = -3.3237, Truncation lag parameter = 4, p-value =
## 0.07095PP.test(mktcap$MKTCAP)##
## Phillips-Perron Unit Root Test
##
## data: mktcap$MKTCAP
## Dickey-Fuller = 0.32774, Truncation lag parameter = 4, p-value = 0.99PP.test(trnvr$TRNR)##
## Phillips-Perron Unit Root Test
##
## data: trnvr$TRNR
## Dickey-Fuller = -4.7199, Truncation lag parameter = 4, p-value = 0.01PP.test(smv$BESN)##
## Phillips-Perron Unit Root Test
##
## data: smv$BESN
## Dickey-Fuller = -0.29909, Truncation lag parameter = 3, p-value =
## 0.9891- Clearly at there exist unit root in all the variables.
Checking PACF
par(mfrow=c(2,2))
pacf(ts(gdpind.zoo, f=1),lag.max = 32, main="PACF of GDP",
xaxt="n", xlab="Lag (quarter)", ylim=c(-.2,1))
axis(1, at=seq(0, 32, 1))
pacf(ts(mktcap.zoo, f=1),lag.max = 32, main="PACF of Market Capitalisation",
xaxt="n", xlab="Lag (quarter)", ylim=c(-.2,1))
axis(1, at=seq(0, 32, 1))
pacf(ts(trnvr.zoo, f=1),lag.max = 32, main="PACF of Turnover",
xaxt="n", xlab="Lag (quarter)", ylim=c(-.2,1))
axis(1, at=seq(0, 32, 1))
pacf(ts(smv.zoo, f=1),lag.max = 32, main="PACF of Stock Volatility",
xaxt="n", xlab="Lag (quarter)", ylim=c(-.2,1))
axis(1, at=seq(0, 32, 1))
- Clearly all are integrated of order 1 , no need to check for stationarity of difference models
Selection of Number of Lagged term by SC Criteria
model=cbind(gdpind.zoo,mktcap.zoo,trnvr.zoo,smv.zoo)
# Writing to a CSV file
write.zoo(x=model,"model.csv")
# Manually removing the NaN values and inserting new values obtained from other sources using
'fix(model)'## [1] "fix(model)"#Importing the updated data set
df1 <- as.data.frame(read_excel(file.path(directory,"updated_values.xlsx"),sheet = "GDP"))
df2 <- as.data.frame(read_excel(file.path(directory,"updated_values.xlsx"),sheet = "MCAP"))
df3 <- as.data.frame(read_excel(file.path(directory,"updated_values.xlsx"),sheet = "TRNVR"))
df4 <- as.data.frame(read_excel(file.path(directory,"updated_values.xlsx"),sheet = "SMV"))
# Converting to zoo series
upd_gdpind.zoo <- zoo(df1$gdpind.zoo,order.by = df1$Index)
upd_mktcap.zoo <- zoo(df2$mktcap.zoo,order.by = df2$Index)
upd_trnvr.zoo <- zoo(df3$trnvr.zoo,order.by = df3$Index)
upd_smv.zoo <- zoo(df4$smv.zoo,order.by = df4$Index)
# Binding all the data into a composite matrix
var_model <- cbind(upd_gdpind.zoo,upd_mktcap.zoo,upd_trnvr.zoo,upd_smv.zoo)
# Deciding the optimal lag value to be included by Schwartz ("SC(n)") criterion
VARselect(var_model,lag.max = 3,type="const")## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 3 2 3
##
## $criteria
## 1 2 3
## AIC(n) 1.139520e+02 1.130875e+02 1.128240e+02
## HQ(n) 1.141592e+02 1.134605e+02 1.133628e+02
## SC(n) 1.144636e+02 1.140084e+02 1.141542e+02
## FPE(n) 3.082240e+49 1.300346e+49 1.002784e+49-Therefore we conclude that the optimal number of lagged terms to be included is : 2
CAUSALITY TEST , APPLYING TODA_YAMAMOTO TEST
$$ \[\begin{aligned} &{\left[\begin{array}{l} M K T C A P_t \\ T R N V R_t \\ S M V_t \\ G D P I N D_t \end{array}\right]=\beta_0+\beta_1\left[\begin{array}{l} M K T C A P_{t-1} \\ T R N V R_{t-1} \\ S M V_{t-1} \\ G D P I N D_{t-1} \end{array}\right]+} \\ &\beta_2\left[\begin{array}{l} M K T C A P_{t-2} \\ T R N V R_{t-2} \\ S M V_{t-2} \\ G D P I N D_{t-2} \end{array}\right]+...\beta_l\left[\begin{array}{l} M K T C A P_{t-l} \\ T R N V R_{t-l} \\ S M V_{t-l} \\ G D P I N D_{t-l} \end{array}\right]+\eta_t \end{aligned}\]$$
\[ For\ Applying\ Toda\ Yamamoto\ test\ we\ need\ to\ construct\ a\ VAR\ model\ with \ k+d_{\max }\ lagged\ terms \ equaling\ 3 \ lags \]
-where k is the number of lagged terms to be included as decided by the SIC criterion and d is the maximum order of integration.
- We may want to do a test of cointegration. If series are cointegrated, there must be a causality. However, Toda and Yamamoto (1995) noted that one advantage of the TY-method is that you don’t have to test for cointegration and, therefore, a pretest bias can be avoided.
ty <- VAR(var_model,p=3,type = "const")
ty##
## VAR Estimation Results:
## =======================
##
## Estimated coefficients for equation upd_gdpind.zoo:
## ===================================================
## Call:
## upd_gdpind.zoo = upd_gdpind.zoo.l1 + upd_mktcap.zoo.l1 + upd_trnvr.zoo.l1 + upd_smv.zoo.l1 + upd_gdpind.zoo.l2 + upd_mktcap.zoo.l2 + upd_trnvr.zoo.l2 + upd_smv.zoo.l2 + upd_gdpind.zoo.l3 + upd_mktcap.zoo.l3 + upd_trnvr.zoo.l3 + upd_smv.zoo.l3 + const
##
## upd_gdpind.zoo.l1 upd_mktcap.zoo.l1 upd_trnvr.zoo.l1 upd_smv.zoo.l1
## 8.749651e-01 1.348220e+06 -3.723653e+07 -4.464410e+08
## upd_gdpind.zoo.l2 upd_mktcap.zoo.l2 upd_trnvr.zoo.l2 upd_smv.zoo.l2
## -1.700111e-01 -1.023550e+06 1.862855e+07 2.986756e+08
## upd_gdpind.zoo.l3 upd_mktcap.zoo.l3 upd_trnvr.zoo.l3 upd_smv.zoo.l3
## 3.337777e-01 -9.386076e+04 1.043888e+07 2.768280e+07
## const
## 3.346709e+11
##
##
## Estimated coefficients for equation upd_mktcap.zoo:
## ===================================================
## Call:
## upd_mktcap.zoo = upd_gdpind.zoo.l1 + upd_mktcap.zoo.l1 + upd_trnvr.zoo.l1 + upd_smv.zoo.l1 + upd_gdpind.zoo.l2 + upd_mktcap.zoo.l2 + upd_trnvr.zoo.l2 + upd_smv.zoo.l2 + upd_gdpind.zoo.l3 + upd_mktcap.zoo.l3 + upd_trnvr.zoo.l3 + upd_smv.zoo.l3 + const
##
## upd_gdpind.zoo.l1 upd_mktcap.zoo.l1 upd_trnvr.zoo.l1 upd_smv.zoo.l1
## -1.050027e-07 1.103332e+00 2.536259e+01 -1.225968e+02
## upd_gdpind.zoo.l2 upd_mktcap.zoo.l2 upd_trnvr.zoo.l2 upd_smv.zoo.l2
## 1.360143e-08 2.412789e-01 -1.892565e+01 1.478810e+02
## upd_gdpind.zoo.l3 upd_mktcap.zoo.l3 upd_trnvr.zoo.l3 upd_smv.zoo.l3
## 1.528581e-07 -5.298624e-01 4.112104e+01 -2.625051e+01
## const
## -6.441721e+05
##
##
## Estimated coefficients for equation upd_trnvr.zoo:
## ==================================================
## Call:
## upd_trnvr.zoo = upd_gdpind.zoo.l1 + upd_mktcap.zoo.l1 + upd_trnvr.zoo.l1 + upd_smv.zoo.l1 + upd_gdpind.zoo.l2 + upd_mktcap.zoo.l2 + upd_trnvr.zoo.l2 + upd_smv.zoo.l2 + upd_gdpind.zoo.l3 + upd_mktcap.zoo.l3 + upd_trnvr.zoo.l3 + upd_smv.zoo.l3 + const
##
## upd_gdpind.zoo.l1 upd_mktcap.zoo.l1 upd_trnvr.zoo.l1 upd_smv.zoo.l1
## -2.947596e-10 1.053096e-03 2.758066e-01 -3.053980e-01
## upd_gdpind.zoo.l2 upd_mktcap.zoo.l2 upd_trnvr.zoo.l2 upd_smv.zoo.l2
## 1.280713e-10 -1.905492e-03 4.621037e-01 1.256720e+00
## upd_gdpind.zoo.l3 upd_mktcap.zoo.l3 upd_trnvr.zoo.l3 upd_smv.zoo.l3
## 5.139362e-10 -4.454741e-04 5.226267e-02 -3.932300e-01
## const
## -3.846356e+03
##
##
## Estimated coefficients for equation upd_smv.zoo:
## ================================================
## Call:
## upd_smv.zoo = upd_gdpind.zoo.l1 + upd_mktcap.zoo.l1 + upd_trnvr.zoo.l1 + upd_smv.zoo.l1 + upd_gdpind.zoo.l2 + upd_mktcap.zoo.l2 + upd_trnvr.zoo.l2 + upd_smv.zoo.l2 + upd_gdpind.zoo.l3 + upd_mktcap.zoo.l3 + upd_trnvr.zoo.l3 + upd_smv.zoo.l3 + const
##
## upd_gdpind.zoo.l1 upd_mktcap.zoo.l1 upd_trnvr.zoo.l1 upd_smv.zoo.l1
## 8.877145e-11 1.770231e-03 7.214014e-03 4.307005e-01
## upd_gdpind.zoo.l2 upd_mktcap.zoo.l2 upd_trnvr.zoo.l2 upd_smv.zoo.l2
## -2.168320e-11 -9.224220e-04 8.258345e-04 2.701003e-01
## upd_gdpind.zoo.l3 upd_mktcap.zoo.l3 upd_trnvr.zoo.l3 upd_smv.zoo.l3
## 4.555919e-11 -5.013238e-04 4.526356e-02 2.808872e-02
## const
## -3.753348e+02# Real GDP (GDPIND) versus Market Capitalization (MKTCAP)
# GDP does not granger cause MKTCAP
wald.test(b=coef(ty$varresult[[1]]), Sigma=vcov(ty$varresult[[1]]), Terms=c(2,6,10))## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 108.2, df = 3, P(> X2) = 0.0 # MKTCAP does not granger cause GDP
wald.test(b=coef(ty$varresult[[2]]), Sigma=vcov(ty$varresult[[2]]), Terms=c(1,5,9))## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 4.6, df = 3, P(> X2) = 0.2# Real GDP (GDPIND) versus Value Traded (TRNVR)
# GDP does not granger cause Value Traded
wald.test(b=coef(ty$varresult[[1]]), Sigma=vcov(ty$varresult[[1]]), Terms=c(3,7,11))## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 6.7, df = 3, P(> X2) = 0.081 # Value Traded does not granger cause GDP
wald.test(b=coef(ty$varresult[[3]]), Sigma=vcov(ty$varresult[[1]]), Terms=c(1,5,9))## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 1.1e-16, df = 3, P(> X2) = 1.0# Real GDP (GDPIND) versus Stock Market Volatility (SMV)
# GDP does not granger cause Stock Market Volatility
wald.test(b=coef(ty$varresult[[1]]), Sigma=vcov(ty$varresult[[1]]), Terms=c(4,8,12))## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 24.2, df = 3, P(> X2) = 2.3e-05 # Stock Market Volatility does not granger GDP
wald.test(b=coef(ty$varresult[[4]]), Sigma=vcov(ty$varresult[[4]]), Terms=c(1,5,9))## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 4.1, df = 3, P(> X2) = 0.25