Test

test.r

#   ------------------------------------------------------------------
#  |PROGRAM NAME: copual and quantmod ex
#  |DATE: 10/30/16
#  |CREATED BY: MATT BOGARD 
#  |PROJECT FILE:              
#  |----------------------------------------------------------------
#  | PURPOSE:         
#  | 
#  | REFERENCE: Modelling Dependence with Copulas in R
#  | http://datascienceplus.com/modelling-dependence-with-copulas/
#  |------------------------------------------------------------------



library(copula)

# use R quantmod to get returns

library(quantmod) # load quantmod

getSymbols("MON")

#get last 3 years of data
mon_dat <- last(MON, '36 months')

head(mon_dat) # view data

# calculate returns

mon_ret <- as.data.frame(dailyReturn(mon_dat)) # returns by day as an R data frame

head(mon_ret) # view data

dim(mon_ret) # check dimensions

# get returns for adm

getSymbols("ADM") # get stock quote data

#get last 3 years of data
adm_dat <- last(ADM, '36 months')

head(adm_dat) # view data

# calculate returns

adm_ret <- as.data.frame(dailyReturn(adm_dat)) # returns by day as an R data frame
print(adm_ret$daily.returns)
head(adm_ret) # view file

dim(adm_ret) # check dimensions



# Plot and correlation checks
plot(mon_ret$daily.returns,adm_ret$daily.returns,pch='.')
abline(lm(adm_ret$daily.returns~mon_ret$daily.returns),col='red',lwd=1)
cor(mon_ret$daily.returns,adm_ret$daily.returns,method='spearman')

# Fit a t-copula
t.cop <- tCopula(dim=2)

set.seed(500)

m <- pobs(as.matrix(cbind(mon_ret$daily.returns,adm_ret$daily.returns)))
fit <- fitCopula(t.cop,m,method='ml')
coef(fit)

# Set the parameters
rho <- coef(fit)[1]
df <- coef(fit)[2]

# Plot the density
persp(tCopula(dim=2,rho,df=df),dCopula)

# Sample random observation from the copula
u <- rCopula(3965,tCopula(dim=2,rho,df=df))
plot(u[,1],u[,2],pch='.',col='blue')
cor(u,method='spearman')

# Estimate marginals parameters
mon_mu <- mean(mon_ret$daily.returns)
mon_sd <- sd(mon_ret$daily.returns)
adm_mu <- mean(adm_ret$daily.returns)
adm_sd <- sd(adm_ret$daily.returns)


# Build the distribution from the copula
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("norm","norm"),
                    paramMargins=list(list(mean=mon_mu, sd=mon_sd),
                                      list(mean=adm_mu, sd=adm_sd)))
# Sample from the distribution
sim <- rmvdc(copula_dist, 3965)


# Compare observed and simulated values
plot(mon_ret$daily.returns,adm_ret$daily.returns,main='Returns')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)