library(tidyquant)
## 필요한 패키지를 로딩중입니다: lubridate
## 
## 다음의 패키지를 부착합니다: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
## 필요한 패키지를 로딩중입니다: PerformanceAnalytics
## 필요한 패키지를 로딩중입니다: xts
## 필요한 패키지를 로딩중입니다: zoo
## 
## 다음의 패키지를 부착합니다: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 다음의 패키지를 부착합니다: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
## 필요한 패키지를 로딩중입니다: quantmod
## 필요한 패키지를 로딩중입니다: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ggplot2)
start_date <- '2018-01-01'
end_date <- '2023-03-31'

sp500_data <- tq_get('SPY',
                     from = start_date,
                     to = end_date,
                     get='stock.prices',
                     periodicity = 'daily')

head(sp500_data)
## # A tibble: 6 x 8
##   symbol date        open  high   low close   volume adjusted
##   <chr>  <date>     <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 SPY    2018-01-02  268.  269.  267.  269. 86655700     245.
## 2 SPY    2018-01-03  269.  271.  269.  270. 90070400     246.
## 3 SPY    2018-01-04  271.  272.  271.  272. 80636400     248.
## 4 SPY    2018-01-05  273.  274.  272.  273. 83524000     249.
## 5 SPY    2018-01-08  273.  274.  273.  274. 57319200     250.
## 6 SPY    2018-01-09  274.  275.  274.  275. 57254000     250.
sp500_return <- sp500_data %>%
  tq_transmute(select = adjusted,
               mutate_fun = dailyReturn,
               type = 'log',
               col_rename = 'return')

head(sp500_return)
## # A tibble: 6 x 2
##   date        return
##   <date>       <dbl>
## 1 2018-01-02 0      
## 2 2018-01-03 0.00631
## 3 2018-01-04 0.00421
## 4 2018-01-05 0.00664
## 5 2018-01-08 0.00183
## 6 2018-01-09 0.00226
sp500_return %>%
  ggplot(aes(x=date,y=return))+
  geom_line(color='#4373B6')+
  theme_bw()+
  labs(x='Date',y='Return',title='S&P500 Return')+
  scale_y_continuous(labels = scales::percent)

ggplot(sp500_return)+
  geom_histogram(aes(x=return, y=after_stat(density)), binwidth = 0.005, color="steelblue", fill="grey", size=1) +
  stat_function(fun = dnorm, args = list(mean = mean(sp500_return$return, na.rm = T), sd = sd(sp500_return$return, na.rm = T)), size=1)+theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

library(PerformanceAnalytics)
library(timetk)
sp500_return_xts <- 
  sp500_return %>% tk_xts()
## Warning: Non-numeric columns being dropped: date
## Using column `date` for date_var.
head(sp500_return_xts)
##                 return
## 2018-01-02 0.000000000
## 2018-01-03 0.006305173
## 2018-01-04 0.004206151
## 2018-01-05 0.006641893
## 2018-01-08 0.001826951
## 2018-01-09 0.002261030
realizedvol <- rollapply(sp500_return_xts,width=22,FUN=sd.annualized)

realizedvol_df <- realizedvol %>% tk_tbl()
tail(realizedvol_df)
## # A tibble: 6 x 2
##   index      return
##   <date>      <dbl>
## 1 2023-03-23  0.177
## 2 2023-03-24  0.178
## 3 2023-03-27  0.178
## 4 2023-03-28  0.174
## 5 2023-03-29  0.180
## 6 2023-03-30  0.180
realizedvol_df %>%
  ggplot(aes(x=index,y=return))+
  geom_line(color='#4373B6')+
  theme_bw()+
  labs(x='Date',y='Volatility',title='S&P 500 Rolling Volatility')+
  scale_y_continuous(labels = scales::percent)
## Warning: Removed 21 rows containing missing values (`geom_line()`).

library(rugarch)
## 필요한 패키지를 로딩중입니다: parallel
## 
## 다음의 패키지를 부착합니다: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
ug_spec <- ugarchspec(mean.model=list(armaOrder=c(1,0)))


ugfit = ugarchfit(spec = ug_spec, data = sp500_return_xts)
ugfit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001033    0.000220   4.6944 0.000003
## ar1    -0.062182    0.030800  -2.0189 0.043499
## omega   0.000005    0.000003   1.8047 0.071120
## alpha1  0.230201    0.030737   7.4893 0.000000
## beta1   0.756602    0.030548  24.7674 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001033    0.000285  3.62456 0.000289
## ar1    -0.062182    0.030160 -2.06176 0.039231
## omega   0.000005    0.000013  0.39607 0.692052
## alpha1  0.230201    0.055151  4.17404 0.000030
## beta1   0.756602    0.088310  8.56758 0.000000
## 
## LogLikelihood : 4152.005 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2833
## Bayes        -6.2637
## Shibata      -6.2834
## Hannan-Quinn -6.2760
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      2.855 0.09107
## Lag[2*(p+q)+(p+q)-1][2]     2.859 0.04634
## Lag[4*(p+q)+(p+q)-1][5]     3.254 0.37018
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.001571  0.9684
## Lag[2*(p+q)+(p+q)-1][5]  0.861658  0.8901
## Lag[4*(p+q)+(p+q)-1][9]  2.343412  0.8603
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]  0.005759 0.500 2.000  0.9395
## ARCH Lag[5]  1.974655 1.440 1.667  0.4769
## ARCH Lag[7]  2.559148 2.315 1.543  0.6003
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.4624
## Individual Statistics:              
## mu     0.12474
## ar1    0.06729
## omega  0.16800
## alpha1 0.08412
## beta1  0.31534
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.28 1.47 1.88
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias           2.7634 0.005801 ***
## Negative Sign Bias  0.8208 0.411909    
## Positive Sign Bias  0.5391 0.589911    
## Joint Effect       13.8664 0.003093 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     60.64    3.068e-06
## 2    30     70.64    2.480e-05
## 3    40     84.24    3.598e-05
## 4    50    109.62    1.559e-06
## 
## 
## Elapsed time : 0.1903021
ugfit@fit$coef
##            mu           ar1         omega        alpha1         beta1 
##  1.033440e-03 -6.218212e-02  5.002358e-06  2.302012e-01  7.566016e-01
ug_var <- ugfit@fit$var %>%ts()
ug_res2 <- (ugfit@fit$residuals)^2 %>% ts()
unloadNamespace('ggfortify')
unloadNamespace('forecast')
library(forecast)

autoplot(ug_var,colour = 'red')+
  autolayer(ug_res2,color = '#4373B6')+
  ggtitle("Plot of Variance and Residual Squared") +
  theme_bw()

ugfore <- ugarchforecast(ugfit,n.ahead = 10)
ugfore
## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: sGARCH
## Horizon: 10
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=2023-03-30]:
##         Series    Sigma
## T+1  0.0007347 0.009702
## T+2  0.0010520 0.009894
## T+3  0.0010323 0.010079
## T+4  0.0010335 0.010259
## T+5  0.0010334 0.010434
## T+6  0.0010334 0.010604
## T+7  0.0010334 0.010768
## T+8  0.0010334 0.010928
## T+9  0.0010334 0.011084
## T+10 0.0010334 0.011235
#Conditional volatility fore next 10 days
ug_f <- ugfore@forecast$sigmaFor %>%ts()


autoplot(ug_f,colour = '#4373B6')+
  labs(x='Day',title='10 Period Forecast of Conditional Volatility')+
  theme_bw()

ug_var_t <- c(tail(ug_var,20),rep(NA,10)) %>% ts()  # gets the last 20 observations
ug_res2_t <- c(tail(ug_res2,20),rep(NA,10)) %>% ts() # gets the last 20 observations
ug_f <- c(rep(NA,20),(ug_f)^2) %>% ts()
ug_f
## Time Series:
## Start = 1 
## End = 30 
## Frequency = 1 
##  [1]           NA           NA           NA           NA           NA
##  [6]           NA           NA           NA           NA           NA
## [11]           NA           NA           NA           NA           NA
## [16]           NA           NA           NA           NA           NA
## [21] 0.0000941231 0.0000978833 0.0001015939 0.0001052555 0.0001088688
## [26] 0.0001124344 0.0001159529 0.0001194250 0.0001228513 0.0001262323
autoplot(ug_res2_t)+
  autolayer(ug_f,series = 'Forecast')+
  autolayer(ug_var_t,series = 'Actual Variance')+
  scale_colour_manual(name='Series', values=c('#4373B6','red'))+
  theme_bw()+
  labs(x='Period',title='Residual^2, Variance, and Forecast')+
   theme(legend.position=c(0.98, 0.98),
        legend.justification=c('right', 'top'))
## Warning: Removed 20 rows containing missing values (`geom_line()`).
## Warning: Removed 10 rows containing missing values (`geom_line()`).