#Getting the data from Yahoo Finance
brka<-getSymbols('BRK-A',src='yahoo',from='2019-01-01',to='2022-02-28',auto.assign=FALSE)
spy<-getSymbols('SPY',src='yahoo',from='2019-01-01',to='2022-02-28',auto.assign=FALSE)
#Subsetting data for training period for the GARCH(1,1) Model
startdate=as.POSIXct('2019-01-01')
enddate=as.POSIXct('2022-01-01')
train_brka=brka[paste(startdate,enddate,sep="::")]
train_spy=spy[paste(startdate,enddate,sep="::")]
chartSeries(brka)
chartSeries(spy)
# Calculating the returns
retb=CalculateReturns(brka$`BRK-A.Adjusted`)
retb=retb[-c(1),]
rets=CalculateReturns(spy$SPY.Adjusted)
rets=rets[-c(1),]
# Calculating the returns
chart_Series(retb)
chart.Histogram(retb,method=c('add.density','add.normal'),colorset = c('blue','red','black'))
legend("topright",legend = c("return","kernel","normal dist"),fill=c('blue','red','black'))
chart_Series(rets)
# Plotting the histogram for the returns
chart.Histogram(rets,method=c('add.density','add.normal'),colorset = c('blue','red','black'))
legend("topright",legend = c("return","kernel","normal dist"),fill=c('blue','red','black'))
#Checking for Normality of the returns
qqnorm(retb,pch=1,frame=FALSE)
qqline(retb,col='steelblue',lwd=2)
qqnorm(rets,pch=1,frame=FALSE)
qqline(rets,col='steelblue',lwd=2)
#Calculating the log returns
logretb=diff(log(brka$`BRK-A.Adjusted`),lag=1)
logrets=diff(log(spy$SPY.Adjusted),lag=1)
chart.Histogram(logretb,method=c('add.density','add.normal'),colorset = c('blue','red','black'))
legend("topright",legend = c("return","kernel","normal dist"),fill=c('blue','red','black'))
chart.Histogram(logrets,method=c('add.density','add.normal'),colorset = c('blue','red','black'))
legend("topright",legend = c("return","kernel","normal dist"),fill=c('blue','red','black'))
#Defining Model Specification for GARCH(1,1), Ftting the log returns for BRKA and forecasting future volatility
model1=ugarchspec(mean.model = list(armaOrder=c(0,0)),variance.model = list(model="sGARCH",garchOrder=c(1,1)),distribution.model='norm')
logretb=na.omit(logretb)# Clearing the NA Values
model_fit=ugarchfit(data=logretb,spec = model1,out.sample=(length(brka$`BRK-A.Adjusted`)-length(train_brka$`BRK-A.Adjusted`)))
model_fit
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000636 0.000354 1.7960 0.072498
## omega 0.000007 0.000002 4.4668 0.000008
## alpha1 0.140958 0.018519 7.6116 0.000000
## beta1 0.807407 0.019771 40.8372 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000636 0.000356 1.7852 0.074233
## omega 0.000007 0.000003 2.7471 0.006013
## alpha1 0.140958 0.017290 8.1525 0.000000
## beta1 0.807407 0.029038 27.8052 0.000000
##
## LogLikelihood : 2329.82
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.1530
## Bayes -6.1285
## Shibata -6.1530
## Hannan-Quinn -6.1435
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.878 0.17055
## Lag[2*(p+q)+(p+q)-1][2] 3.672 0.09237
## Lag[4*(p+q)+(p+q)-1][5] 5.446 0.12110
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.3418 0.5588
## Lag[2*(p+q)+(p+q)-1][5] 1.0336 0.8520
## Lag[4*(p+q)+(p+q)-1][9] 2.1045 0.8922
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1455 0.500 2.000 0.7029
## ARCH Lag[5] 1.4475 1.440 1.667 0.6062
## ARCH Lag[7] 1.8175 2.315 1.543 0.7559
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 4.7729
## Individual Statistics:
## mu 0.04653
## omega 1.25881
## alpha1 0.20719
## beta1 0.11140
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.608 0.1083
## Negative Sign Bias 1.466 0.1430
## Positive Sign Bias 1.008 0.3139
## Joint Effect 3.544 0.3151
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 15.11 0.7155
## 2 30 32.10 0.3157
## 3 40 43.47 0.2867
## 4 50 54.05 0.2875
##
##
## Elapsed time : 0.108212
forecast_brka<-ugarchforecast(model_fit,n.ahead=38)
forecast_brka
##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: sGARCH
## Horizon: 38
## Roll Steps: 0
## Out of Sample: 38
##
## 0-roll forecast [T0=2021-12-31]:
## Series Sigma
## T+1 0.0006364 0.009317
## T+2 0.0006364 0.009453
## T+3 0.0006364 0.009580
## T+4 0.0006364 0.009699
## T+5 0.0006364 0.009810
## T+6 0.0006364 0.009915
## T+7 0.0006364 0.010013
## T+8 0.0006364 0.010105
## T+9 0.0006364 0.010192
## T+10 0.0006364 0.010273
## T+11 0.0006364 0.010350
## T+12 0.0006364 0.010422
## T+13 0.0006364 0.010490
## T+14 0.0006364 0.010554
## T+15 0.0006364 0.010615
## T+16 0.0006364 0.010672
## T+17 0.0006364 0.010726
## T+18 0.0006364 0.010776
## T+19 0.0006364 0.010824
## T+20 0.0006364 0.010869
## T+21 0.0006364 0.010912
## T+22 0.0006364 0.010953
## T+23 0.0006364 0.010991
## T+24 0.0006364 0.011027
## T+25 0.0006364 0.011061
## T+26 0.0006364 0.011093
## T+27 0.0006364 0.011124
## T+28 0.0006364 0.011152
## T+29 0.0006364 0.011180
## T+30 0.0006364 0.011205
## T+31 0.0006364 0.011230
## T+32 0.0006364 0.011253
## T+33 0.0006364 0.011275
## T+34 0.0006364 0.011296
## T+35 0.0006364 0.011315
## T+36 0.0006364 0.011334
## T+37 0.0006364 0.011351
## T+38 0.0006364 0.011368
variance_rate_brka<-coef(model_fit)[2]/(1-coef(model_fit)[3]-coef(model_fit)[4])
LT_vol_brka=sqrt(coef(model_fit)[2]/(1-coef(model_fit)[3]-coef(model_fit)[4]))
#Defining Model Specification for GARCH(1,1), Ftting the log returns for SPY and forecasting future volatility
model2=ugarchspec(mean.model=list(armaOrder=c(0,0)),variance.model=list(model='sGARCH',garchOrder=c(1,1)),distribution.model = 'norm')
logrets=na.omit(logrets)
model2_fit=ugarchfit(data=logrets,spec=model2,out.sample = (length(spy$SPY.Adjusted)-length(train_spy$SPY.Adjusted)))
model2_fit
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001272 0.000268 4.7463 0.000002
## omega 0.000007 0.000002 3.8050 0.000142
## alpha1 0.285623 0.042327 6.7480 0.000000
## beta1 0.682136 0.022975 29.6909 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001272 0.000247 5.1493 0.000000
## omega 0.000007 0.000003 1.9791 0.047800
## alpha1 0.285623 0.056981 5.0126 0.000001
## beta1 0.682136 0.067770 10.0654 0.000000
##
## LogLikelihood : 2461.985
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.5026
## Bayes -6.4781
## Shibata -6.5027
## Hannan-Quinn -6.4932
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.246 0.2643
## Lag[2*(p+q)+(p+q)-1][2] 1.694 0.3190
## Lag[4*(p+q)+(p+q)-1][5] 2.745 0.4555
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.4782 0.4892
## Lag[2*(p+q)+(p+q)-1][5] 1.4406 0.7545
## Lag[4*(p+q)+(p+q)-1][9] 2.4803 0.8406
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.01209 0.500 2.000 0.9124
## ARCH Lag[5] 1.39278 1.440 1.667 0.6209
## ARCH Lag[7] 1.64297 2.315 1.543 0.7923
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.5416
## Individual Statistics:
## mu 0.07789
## omega 0.32591
## alpha1 0.15593
## beta1 0.13604
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.7064 0.08835 *
## Negative Sign Bias 0.4286 0.66832
## Positive Sign Bias 0.5488 0.58328
## Joint Effect 5.4766 0.14004
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 40.67 0.0026725
## 2 30 61.14 0.0004443
## 3 40 68.44 0.0024629
## 4 50 85.53 0.0009574
##
##
## Elapsed time : 0.09510994
forecast_spy<-ugarchforecast(model2_fit,n.ahead=38)
forecast_spy
##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: sGARCH
## Horizon: 38
## Roll Steps: 0
## Out of Sample: 38
##
## 0-roll forecast [T0=2021-12-31]:
## Series Sigma
## T+1 0.001272 0.007256
## T+2 0.001272 0.007588
## T+3 0.001272 0.007897
## T+4 0.001272 0.008184
## T+5 0.001272 0.008453
## T+6 0.001272 0.008705
## T+7 0.001272 0.008942
## T+8 0.001272 0.009166
## T+9 0.001272 0.009377
## T+10 0.001272 0.009578
## T+11 0.001272 0.009768
## T+12 0.001272 0.009948
## T+13 0.001272 0.010119
## T+14 0.001272 0.010282
## T+15 0.001272 0.010438
## T+16 0.001272 0.010586
## T+17 0.001272 0.010728
## T+18 0.001272 0.010863
## T+19 0.001272 0.010992
## T+20 0.001272 0.011116
## T+21 0.001272 0.011234
## T+22 0.001272 0.011348
## T+23 0.001272 0.011456
## T+24 0.001272 0.011561
## T+25 0.001272 0.011661
## T+26 0.001272 0.011757
## T+27 0.001272 0.011849
## T+28 0.001272 0.011937
## T+29 0.001272 0.012022
## T+30 0.001272 0.012104
## T+31 0.001272 0.012182
## T+32 0.001272 0.012258
## T+33 0.001272 0.012330
## T+34 0.001272 0.012400
## T+35 0.001272 0.012468
## T+36 0.001272 0.012532
## T+37 0.001272 0.012595
## T+38 0.001272 0.012655
variance_rate_spy<-coef(model_fit)[2]/(1-coef(model_fit)[3]-coef(model_fit)[4])
LT_vol_spy=sqrt(coef(model2_fit)[2]/(1-coef(model2_fit)[3]-coef(model2_fit)[4]))
# VaR Estimation based on returns
retb<-(train_brka$`BRK-A.Adjusted`-lag.xts(train_brka$`BRK-A.Adjusted`))/lag.xts(train_brka$`BRK-A.Adjusted`)
rets<-(train_spy$SPY.Adjusted-lag.xts(train_spy$SPY.Adjusted))/lag.xts(train_spy$SPY.Adjusted)
retb=na.omit(retb)
rets=na.omit(rets)
corr<-cor(retb,rets)
#Daily Positional Losses historically
portfolio_size=1000000
daily_pos_b=retb*portfolio_size
daily_pnl_b=daily_pos_b
daily_pnl_b=na.omit(daily_pnl_b)
VaR_brka=quantile(daily_pnl_b,0.01)
ES_brka<-mean(sort(as.numeric(daily_pnl_b))[1:7])
daily_pos_s=rets*portfolio_size
daily_pnl_s=daily_pos_s
daily_pnl_s=na.omit(daily_pnl_s)
VaR_spy=quantile(daily_pnl_s,0.01)
ES_spy<-mean(sort(as.numeric(daily_pnl_s))[1:7])
daily_pos_port=daily_pos_s+daily_pos_b
daily_pnl_port=na.omit(daily_pos_port)
VaR_port=quantile(daily_pnl_port,0.01)
ES_port=mean(sort(as.numeric(daily_pnl_port))[1:7])
diversification_benefit=VaR_port-VaR_spy-VaR_brka
#Plotting the positional losses
temp1<-sort(as.numeric(daily_pnl_b))
temp2<-sort(as.numeric(daily_pnl_s))
temp3<-sort(as.numeric(daily_pnl_s+daily_pnl_b))
pnl<-data.frame(cbind(temp1,temp2,temp3))
colnames(pnl)<-c('PnL_BRKA','PnL_SPY','PnL_Portfolio')
g<-ggplot(pnl,mapping=aes(x=PnL_BRKA))
g+geom_histogram(color='darkblue',fill='lightblue')+geom_vline(data=pnl,aes(xintercept=quantile(PnL_BRKA,0.01),col='VaR'))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
g2<-ggplot(pnl,mapping=aes(x=PnL_SPY))
g2+geom_histogram(color='darkblue',fill='lightblue')+geom_vline(data=pnl,aes(xintercept=quantile(PnL_SPY,0.01),col='VaR'))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
g3<-ggplot(pnl,mapping=aes(x=PnL_Portfolio))
g3+geom_histogram(color='darkblue',fill='lightblue')+geom_vline(data=pnl,aes(xintercept=quantile(PnL_Portfolio,0.01),col='VaR'))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#VaR backtesting:
#Subset the dates to get the data
backtest_start=as.POSIXct('2021-02-28')
backtest_end=as.POSIXct('2022-02-28')
backtest_brka=brka[paste(backtest_start,backtest_end,sep="::")]
backtest_spy=spy[paste(backtest_start,backtest_end,sep="::")]
#Calculating the daily returns for the period
retb_bt<-CalculateReturns(backtest_brka$`BRK-A.Adjusted`)
rets_bt<-CalculateReturns(backtest_spy$SPY.Adjusted)
daily_pos_brka_bt<-rets_bt*portfolio_size
daily_pos_spy_bt<-rets_bt*portfolio_size
daily_pos_brka_bt=na.omit(daily_pos_brka_bt)
daily_pos_spy_bt=na.omit(daily_pos_spy_bt)
#Number of days the loss crosses VaR
count_spy=sum(daily_pos_spy_bt<VaR_spy)
count_spy
## [1] 0
count_brka=sum(daily_pos_brka_bt<VaR_brka)
count_brka
## [1] 0
daily_pos_port_bt=daily_pos_brka_bt+daily_pos_spy_bt
count_port=sum(daily_pos_port_bt<VaR_port)
count_port
## [1] 0
#VaR BASED ON NORMAL DISTRIBUTION
sigb=sd(retb)
sigs=sd(rets)
portsig=(sigb^2/4+sigs^2/4+2*sigs*sigb*corr)^0.5
VaR_b_N=-qnorm(0.99)*sigb*10^6
VaR_s_N=-qnorm(0.99)*sigs*10^6
VaR_p_N=-qnorm(0.99)*portsig*2*10^6
VaR_b_N
## [1] -33150.44
VaR_s_N
## [1] -32121.53
VaR_p_N
## SPY.Adjusted
## BRK-A.Adjusted -94621
#ECONOMIC CAPITAL
RAROC_spy<-as.numeric(rets_bt[paste(as.POSIXct('2022-01-03'))])*10^6/VaR_spy
RAROC_brka<-as.numeric(retb_bt[paste(as.POSIXct('2022-01-03'))])*10^6/VaR_brka
RAROC_port<-(as.numeric(rets_bt[paste(as.POSIXct('2022-01-03'))])+as.numeric(retb_bt[paste(as.POSIXct('2022-01-03'))]))*10^6/VaR_port
RAROC_brka
## 1%
## -0.2096993
RAROC_spy
## 1%
## -0.1318538
RAROC_port
## 1%
## -0.1679277
#export the data for the report
#write.csv(as.data.frame(brka),"brka.csv")
#write.csv(as.data.frame(spy),"spy.csv")