library(ggplot2,forecast)
library(astsa)
library(zoo,lmtest)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(fUnitRoots)
## Loading required package: timeDate
## Loading required package: timeSeries
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
## Loading required package: fBasics
## 
## Attaching package: 'fBasics'
## The following object is masked from 'package:astsa':
## 
##     nyse
library(FitARMA)
## Loading required package: FitAR
## Loading required package: lattice
## Loading required package: leaps
## Loading required package: ltsa
## Loading required package: bestglm
library(strucchange)
## Loading required package: sandwich
library(reshape)
library(Rmisc)
## Loading required package: plyr
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:reshape':
## 
##     rename, round_any
library(fBasics)
library(tsoutliers)
library(TSA)
## Loading required package: locfit
## locfit 1.5-9.1    2013-03-22
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-20. For overview type 'help("mgcv-package")'.
## Loading required package: tseries
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:timeDate':
## 
##     kurtosis, skewness
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(dygraphs)
library(quantmod)
## Loading required package: xts
## Loading required package: TTR
## 
## Attaching package: 'TTR'
## The following object is masked from 'package:fBasics':
## 
##     volatility
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:plyr':
## 
##     here
## The following object is masked from 'package:reshape':
## 
##     stamp
## The following object is masked from 'package:base':
## 
##     date

{r warning=FALSE} # dji <- read.csv("C:\\Users\\r631758\\Desktop\\r631758\\R codes\\TimeSeries\\data\\^DJI.csv", header=TRUE) # SP <- read.csv("C:\\Users\\r631758\\Desktop\\r631758\\R codes\\TimeSeries\\data\\^GSPC.csv", header=TRUE) # IXIC <- read.csv("C:\\Users\\r631758\\Desktop\\r631758\\R codes\\TimeSeries\\data\\^IXIC.csv", header=TRUE) #

tickers <- c("^DJI", "^GSPC", "^IXIC")

getSymbols(tickers, from=as.Date("1950-01-15"), to=as.Date("2018-01-16"))
## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.
## 
## WARNING: There have been significant changes to Yahoo Finance data.
## Please see the Warning section of '?getSymbols.yahoo' for details.
## 
## This message is shown once per session and may be disabled by setting
## options("getSymbols.yahoo.warning"=FALSE).
## [1] "DJI"  "GSPC" "IXIC"
symbol <- c("DJI", "GSPC", "IXIC")
closePrices <- do.call(merge, lapply(symbol, function(x) Cl(get(x))))
dateWindowDJI <- c("1985-01-29", "2018-01-16")
dateWindowGSPS <- c("1950-01-16", "2018-01-16")
dateWindowIXIC <- c("1971-02-05", "2018-01-16")

dygraph(closePrices$DJI.Close["1985-01-29/2018"], main = "DJI", group = "stock") %>%
    dyRangeSelector(dateWindow = dateWindowDJI)%>% dyEvent("1999-06-30", "crisis1", labelLoc = "bottom")
dygraph(closePrices$GSPC.Close, main = "GSPC", group = "stock") %>%
    dyRangeSelector(dateWindow = dateWindowGSPS)
dygraph(closePrices$IXIC.Close, main = "IXIC", group = "stock") %>%
    dyRangeSelector(dateWindow = dateWindowIXIC)
# dygraph(closePrices, main = "Percent", group = "stock") %>%
#   dyRebase(percent = TRUE) %>%
#   dyRangeSelector(dateWindow = dateWindow)
# 
# dygraph(closePrices, main = "None", group = "stock") %>%
#   dyRangeSelector(dateWindow = dateWindow)
DJI<-closePrices$DJI.Close["1985-01-29/2018-01-16"]
tt<-1:length(DJI)
fit<-ts(loess(DJI~tt, span=0.2)$fitted, frequency = 1)
#assign date to the fit data
d<-index(DJI)
fit1<-xts(fit,as.Date(d))
DJI<-closePrices$DJI.Close["1985-01-29/2018"]
test<-cbind(DJI, fit1)


dygraph(test,main = "Value", group = "stock") 
plot(density(DJI))

run Dickey-Fuller test

urdftest_lag=floor(12*(length(DJI)/100)^0.25)
urdfTest(DJI,lags = urdftest_lag, type = c("ct"), doplot = FALSE)
## 
## Title:
##  Augmented Dickey-Fuller Unit Root Test
## 
## Test Results:
## 
## 
## Description:
##  Tue Jan 16 14:43:41 2018 by user: r631758
urkpssTest(DJI, type = c("tau"), lags = c("long"),  doplot = FALSE)
## 
## Title:
##  KPSS Unit Root Test
## 
## Test Results:
##   NA
## 
## Description:
##  Tue Jan 16 14:43:41 2018 by user: r631758
summary(lm(DJI ~ 1))
## 
## Call:
## lm(formula = DJI ~ 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7642.0 -5437.8   668.7  3099.2 16919.2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8884.02      58.63   151.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5345 on 8308 degrees of freedom
lev_fit <- lm(DJI ~ 1)
summary(lev_fit)
## 
## Call:
## lm(formula = DJI ~ 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7642.0 -5437.8   668.7  3099.2 16919.2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8884.02      58.63   151.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5345 on 8308 degrees of freedom

search structural changes

DJI_window<-DJI["2008/2009"]

#DJI_window<-DJI["2008-01-20/2008-03-23"]
d<-index(DJI_window)
DJI_window1<-ts(DJI_window$DJI.Close)
dateMap<-cbind.data.frame(index(DJI_window1), as.Date(index(DJI_window)))


DJI_brk <- breakpoints(DJI_window1~ 1, h = 0.1)
DJI_brk$breakpoints<-dateMap$`as.Date(index(DJI_window))`[c(DJI_brk$breakpoints)]
summary(DJI_brk)
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = DJI_window1 ~ 1, h = 0.1)
## 
## Breakpoints at observation number:
##                                           
## m = 1              191                    
## m = 2              192             393    
## m = 3      118     193             393    
## m = 4      118     192 261         392    
## m = 5      118     192 264     333 399    
## m = 6      118     192 264     333 392 447
## m = 7   61 118     192 264     333 392 447
## m = 8   61 118     192 242 292 342 392 447
## m = 9   50 100 150 200 255 305 355 405 455
## 
## Corresponding to breakdates:
##                                           
## m = 1              191                    
## m = 2              192             393    
## m = 3      118     193             393    
## m = 4      118     192 261         392    
## m = 5      118     192 264     333 399    
## m = 6      118     192 264     333 392 447
## m = 7   61 118     192 264     333 392 447
## m = 8   61 118     192 242 292 342 392 447
## m = 9   50 100 150 200 255 305 355 405 455
## 
## Fit:
##                                                                          
## m   0         1         2         3         4         5         6        
## RSS 1.562e+09 3.502e+08 1.730e+08 1.089e+08 8.943e+07 6.779e+07 5.337e+07
## BIC 8.993e+03 8.250e+03 7.906e+03 7.685e+03 7.598e+03 7.471e+03 7.362e+03
##                                  
## m   7         8         9        
## RSS 5.178e+07 6.243e+07 1.023e+08
## BIC 7.359e+03 7.466e+03 7.728e+03
P1<-plot(DJI_brk)

DJI_brk$breakpoints
## [1] "2008-03-31" "2008-06-19" "2008-10-03" "2009-01-16" "2009-04-28"
## [6] "2009-07-22" "2009-10-08"
coef(DJI_brk, breaks = length(DJI_brk$breakpoints))
##           (Intercept)
## 1 - 61      12386.385
## 62 - 118    12618.366
## 119 - 192   11320.644
## 193 - 264    8707.421
## 265 - 333    7680.782
## 334 - 392    8484.017
## 393 - 447    9469.375
## 448 - 505   10226.000
DJI_brk1<-xts(fitted(DJI_brk,breaks=length(DJI_brk$breakpoints)),as.Date(d))



Breakplot<-cbind(DJI_window,DJI_brk1)
dygraph(Breakplot) %>%dyRangeSelector(dateWindow = c(DJI_brk$breakpoints[1], DJI_brk$breakpoints[length(DJI_brk$breakpoints)]))%>%dyAxis("x",drawGrid = FALSE)%>%dyAxis("y",drawGrid = FALSE)
CI<-confint(DJI_brk,breaks=length(DJI_brk$breakpoints))

LOWCI<-dateMap$`as.Date(index(DJI_window))`[match(CI$confint[,1], dateMap$`index(DJI_window1)`)]
middle<-dateMap$`as.Date(index(DJI_window))`[match(CI$confint[,2], dateMap$`index(DJI_window1)`)]           
HCI<-dateMap$`as.Date(index(DJI_window))`[match(CI$confint[,3], dateMap$`index(DJI_window1)`)]  

CIDate<-cbind.data.frame(LOWCI,middle, HCI)

CIDate
plot(DJI_window1)
lines(fitted(DJI_brk,breaks=length(DJI_brk$breakpoints)),col=4)
lines(confint(DJI_brk,breaks=length(DJI_brk$breakpoints)))

trend structual changes

DJI_window<-DJI["2008/2009"]

#DJI_window<-DJI["2008-01-20/2008-03-23"]
d<-index(DJI_window)
DJI_window1<-ts(DJI_window$DJI.Close)
dateMap<-cbind.data.frame(index(DJI_window1), as.Date(index(DJI_window)))

tt<-1:length(DJI_window)
trend_fit <- lm(DJI_window1~ tt)
summary(trend_fit)
## 
## Call:
## lm(formula = DJI_window1 ~ tt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3170.5 -1259.8   408.8   986.5  2431.8 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12033.1735   120.1588  100.14   <2e-16 ***
## tt             -7.7707     0.4115  -18.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1348 on 503 degrees of freedom
## Multiple R-squared:  0.4148, Adjusted R-squared:  0.4137 
## F-statistic: 356.6 on 1 and 503 DF,  p-value: < 2.2e-16
plot(DJI_window1)
lines(ts(fitted(trend_fit),start=1, frequency=1), col = 4)

DJI_brk<-breakpoints(DJI_window1~tt, h=0.1)


DJI_brk$breakpoints<-dateMap$`as.Date(index(DJI_window))`[c(DJI_brk$breakpoints)]
summary(DJI_brk)
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = DJI_window1 ~ tt, h = 0.1)
## 
## Breakpoints at observation number:
##                                           
## m = 1                      283            
## m = 2              192     283            
## m = 3      108     192     283            
## m = 4      108     192 251 301            
## m = 5      83  136 193 251 301            
## m = 6      83  136 193 251 301 365        
## m = 7      83  136 193 251 301 370 420    
## m = 8      83  136 192 242 292 342 392 447
## m = 9   50 100 150 200 250 300 350 400 450
## 
## Corresponding to breakdates:
##                                           
## m = 1                      283            
## m = 2              192     283            
## m = 3      108     192     283            
## m = 4      108     192 251 301            
## m = 5      83  136 193 251 301            
## m = 6      83  136 193 251 301 365        
## m = 7      83  136 193 251 301 370 420    
## m = 8      83  136 192 242 292 342 392 447
## m = 9   50 100 150 200 250 300 350 400 450
## 
## Fit:
##                                                                          
## m   0         1         2         3         4         5         6        
## RSS 914155475 159094390  57697906  39215807  34613106  30109241  28556371
## BIC      8728      7864      7370      7194      7150      7098      7090
##                                  
## m   7         8         9        
## RSS  26580223  27258432  35362099
## BIC      7072      7104      7254
P1<-plot(DJI_brk)

DJI_brk$breakpoints
## [1] "2008-04-30" "2008-07-16" "2008-10-06" "2008-12-29" "2009-03-12"
## [6] "2009-06-19" "2009-08-31"
coef(DJI_brk, breaks = length(DJI_brk$breakpoints))
##           (Intercept)          tt
## 1 - 83      12419.817   0.9104084
## 84 - 136    16720.502 -41.2861144
## 137 - 193   13952.436 -16.1931378
## 194 - 251   10876.213  -9.7649505
## 252 - 301   20200.364 -44.6560227
## 302 - 370    1659.664  19.3884483
## 371 - 420   -2926.795  29.9883384
## 421 - 505    3858.777  13.3473942
DJI_brk1<-xts(fitted(DJI_brk,breaks=length(DJI_brk$breakpoints)),as.Date(d))



Breakplot<-cbind(DJI_window,DJI_brk1)
dygraph(Breakplot) %>%dyRangeSelector(dateWindow = c(DJI_brk$breakpoints[1], DJI_brk$breakpoints[length(DJI_brk$breakpoints)]))%>%dyAxis("x",drawGrid = FALSE)%>%dyAxis("y",drawGrid = FALSE)
CI<-confint(DJI_brk,breaks=length(DJI_brk$breakpoints))

LOWCI<-dateMap$`as.Date(index(DJI_window))`[match(CI$confint[,1], dateMap$`index(DJI_window1)`)]
middle<-dateMap$`as.Date(index(DJI_window))`[match(CI$confint[,2], dateMap$`index(DJI_window1)`)]           
HCI<-dateMap$`as.Date(index(DJI_window))`[match(CI$confint[,3], dateMap$`index(DJI_window1)`)]  

CIDate<-cbind.data.frame(LOWCI,middle, HCI)

CIDate
plot(DJI_window1)
lines(fitted(DJI_brk,breaks=length(DJI_brk$breakpoints)),col=4)
lines(confint(DJI_brk,breaks=length(DJI_brk$breakpoints)))