rm(list = ls());gc()
library(data.table)
library(copula)
library(tidyverse)
library(quantmod)
library(PerformanceAnalytics)

# 抓股價報酬率資料

getStock<-function(code){
    stock<-getSymbols(paste0(code,'.TW'),from='2010-1-1',to='2019-12-31',auto.assign = F)
    data<-stock %>% 
    as.data.frame(row.names = NULL)%>%
    rename_all(~c('Open','High','Low','Close','Volume','Adjusted')) %>%
    mutate(date = index(stock)) %>% transmute(date,rt = Adjusted / shift(Adjusted) - 1) %>% 
    na.omit() %>%
    rename_at(vars(rt),~code)
    return(data)
}
s2330<-getStock('2330')
s2317<-getStock('2317')

# 抓匯率資料

TWD <- getSymbols('TWD=X',from = '2010-1-1',to = '2019-12-31', auto.assign = F)
TWDdata <- TWD%>% 
  as.data.frame(row.names = NULL) %>%
  rename_all(~c('Open','High','Low','Close','Volume','Adjusted')) %>%
  mutate(date = index(TWD)) %>% 
  select(date, Adjusted) %>% 
  rename_at(vars(Adjusted), ~'TWD')

CNY <- getSymbols('CNY=X',from = '2010-1-1',to = '2019-12-31', auto.assign = F)
CNYdata <- CNY %>% 
  as.data.frame(row.names = NULL) %>%
  rename_all(~c('Open','High','Low','Close','Volume','Adjusted')) %>%
  mutate(date = index(TWD)) %>% 
  select(date, Adjusted) %>% 
  rename_at(vars(Adjusted), ~'CNY')

USDTWD <- TWDdata %>%
    transmute(date, USDTWD = TWD / shift(TWD) - 1) %>% na.omit()
CNYTWD <- TWDdata %>%
    left_join(CNYdata, by = 'date') %>%
    mutate(cross = TWD/CNY) %>%
    transmute(date, CNYTWD = cross / shift(cross) - 1) %>% na.omit()

# 合併

Data<-s2330 %>% inner_join(s2317,by='date')%>%
  inner_join(USDTWD,by='date')%>%
  inner_join(CNYTWD,by='date')%>% 
  mutate(date=NULL)

# copula

normal.cop <- normalCopula(dim = 4)
fit.cop <- fitCopula(normal.cop, pobs(as.matrix(Data)), method = 'ml')
# create coef
rho <- coef(fit.cop)
# print(rho)

# 模擬

set.seed(100)
u1 <- rCopula(10000,normalCopula(rho, dim = 4))

# 替換成原本資料

u2=u1
for (a in 1:4) {
  for (b in 1:nrow(u2)) {
    u2[b,a]=ceiling(u2[b,a]*nrow(Data))
  }
}
u2[,1]<-Data[u2[,1],1]
u2[,2]<-Data[u2[,2],2]
u2[,3]<-Data[u2[,3],3]
u2[,4]<-Data[u2[,4],4]

# 計算VaR

port<-ts(0.2*x.sample[,1]+0.2*x.sample[,2]+0.4*x.sample[,3]+0.2*x.sample[,4])
VaR(port,0.95)
quantile(port, 0.5)-quantile(port, 0.05)

VaR -0.009195181
quantile 0.009597874