# ! Set your working directory (folder containing the subfolders R_allocation, R_interpretation, data, weights, etc)
setwd("c:/Documents and Settings/Administrator/Desktop/risk budget programs")
# Programs:
source("R_Allocation/allocation_functions_monthly.R"); source("R_Allocation/estimators.R");
library(zoo); library(fGarch); library("PerformanceAnalytics");
memory.limit(2048)
# Choose data
datacase = "equitybondscommodity" # "equitybondscommodity" or "equitybondscommodity"
# Define optimization criteria
criteria = c( "EW" )
# Yearly or quarterly rebalancing
yearly = F;
# Length estimation period
estyears = 4
source("R_Allocation/allocation_functions_monthly.R"); source("R_Allocation/estimators.R");
library(zoo); library(fGarch); library("PerformanceAnalytics");
memory.limit(2048)
#############################################################################
# Load the data
#############################################################################
firstyear = 1976 ; firstquarter = 1; lastyear = 2009; lastquarter = 2;
data = read.table( file= paste("data/",datacase,"/data.txt",sep="") ,header=T)
date = as.Date(data[,1],format="%Y-%m-%d")
data = zoo( data[,2:ncol(data)] , order.by = date )
# "Bond" "SP500" "EAFE" "SPGSCI" "TBill" "Inflation"
cAssets = ncol(data)-2; # number of risky assets
# The loaded data has monthly frequency
# Estimation uses real returns
monthlyR = data[,(1:(cAssets+1))]-data[,cAssets+2]
monthlyR.cash = monthlyR[,cAssets+1]
R = monthlyR[,1:cAssets]
plot(R)
head(R);tail(R)
from=to=c();
if(yearly){
for( year in firstyear:(lastyear-estyears+1) ){
from = c( from , paste( year,"-01-01",sep="" ) );
to = c( to , paste( (year+estyears-1),"-12-31",sep="" ) ) }
nsamples = length(from);
names.input = paste( rep("y_",nsamples) , seq( (firstyear+estyears),lastyear+1,1) , sep="" );
}else{
for( year in firstyear:(lastyear-estyears) ){
from = c( from , paste( year,"-01-01",sep="" ) , paste( year,"-04-01",sep="" ) , paste( year,"-07-01",sep="" ) , paste( year,"-10-01",sep="" ) );
to = c( to , paste( (year+estyears-1),"-12-31",sep="" ), paste( (year+estyears),"-03-31",sep="" ), paste( (year+estyears),"-06-30",sep="" ),
paste( (year+estyears),"-09-30",sep="" ) ) }
from = c( from , paste( year+1,"-01-01",sep="" ) ); to = c( to , paste( (year+estyears),"-12-31",sep="" ) )
nsamples = length(from);
#names of quarters for which the forecast is made:
names.input = paste( c("Q1y_","Q2y_","Q3y_","Q4y_") , rep(seq( (firstyear+estyears),lastyear,1),each=4) , sep="" );
names.input = c( names.input , paste( "Q1y_" , lastyear+1, sep="" ) );
# this is assuming estimation window runs from firstyearQ1 to lastyearQ4
from = from[firstquarter:(length(from)-4+lastquarter)]; to = to[firstquarter:(length(to)-4+lastquarter)]
names.input = names.input[firstquarter:(length(names.input)-4+lastquarter)]
}
cPeriods = length(to) ; cCriteria = length(criteria) ;
percriskcontribcriterion = "mES";
RCout = matrix( rep(0,cPeriods*cAssets) , ncol= cAssets );
#wper = matrix(rep(1/cAssets, cAssets),ncol=1)
W = read.csv( file = paste("C:\\Documents and Settings\\Administrator\\Desktop\\risk budget programs\\weights\\mES\\equitybondscommodity\\quarterly\\unconstrained\\", "mES.4yr-InfInf" , ".csv" , sep="") );
RCout = matrix( rep(0, cCriteria*cPeriods*(cAssets)) , ncol= (cAssets) ); # add possibility of cash if risk budget constraints are infeasible.
# downside risk
alpha = 0.05;
# Estimation of the return mean vector, covariance, coskewness and cokurtosis matrix
for( per in c(1:cPeriods) ){
print("-----------New estimation period ends on observation------------------")
print( paste(to[per],"out of total number of obs equal to", max(to) ));
print("----------------------------------------------------------------")
matrix(as.vector(wper[1:4]),ncol=1)
# Estimate GARCH model with data from inception
inception.R = window(R, start = as.Date(from[1]) , end = as.Date(to[per]) );
# Estimate comoments of innovations with rolling estimation windows
in.sample.R = window(R, start = as.Date(from[per]) , end = as.Date(to[per]) );
in.sample.R = checkData(in.sample.R, method="matrix");
# Estimation of mean return
M = c();
library(TTR)
Tmean = 47 # monthly returns: 4 year exponentially weighted moving average
for( i in 1:cAssets ){
# M = cbind( M , as.vector( rollmean(x=inception.R[,i],k=Tmean, align="right",na.pad=T) ) ) #2/(n+1)
M = cbind( M , as.vector( EMA(x=inception.R[,i],n=Tmean) ) ) #2/(n+1)
# M = cbind( M , as.vector(cumsum())/c(1:length(inception.R[,i])) )
}
M = zoo( M , order.by=time(inception.R) )
# Center returns (shift by one observations since M[t,] is rolling mean t-Tmean+1,...,t; otherwise lookahead bias)
inception.R.cent = inception.R;
ZZ = matrix( rep(as.vector( apply( inception.R[1:Tmean, ] , 2 , 'mean' )),Tmean),byrow=T,nrow=Tmean);
inception.R.cent[1:Tmean,] = inception.R[1:Tmean, ] - ZZ;
if( nrow(inception.R)>(Tmean+1) ){
A = M[Tmean:(nrow(inception.R)-1),];
A = zoo( A , order.by = time(inception.R[(Tmean+1):nrow(inception.R), ])) ; #shift dates; otherwise zoo poses problem
inception.R.cent[(Tmean+1):nrow(inception.R), ] = inception.R[(Tmean+1):nrow(inception.R), ] - A}
# Garch estimation
S = c();
for( i in 1:cAssets ){
gout = garchFit(formula ~ garch(1,1), data = inception.R.cent[,i],include.mean = F, cond.dist="QMLE", trace = FALSE )
if( as.vector(gout@fit$coef["alpha1"]) < 0.01 ){
sigmat = rep( sd( as.vector(inception.R.cent[,i])), length(inception.R.cent[,i]) );
}else{
sigmat = gout@sigma.t
}
S = cbind( S , sigmat)
}
S = zoo( S , order.by=time(inception.R.cent) )
#plot(S, main="Check visually garch volatility estimates \n If estimated alpha < 0.01, the standard deviation is used. " )
# Estimate correlation, coskewness and cokurtosis matrix locally using cleaned innovation series in three year estimation window
selectU = window(inception.R.cent, start = as.Date(from[per]) , end = as.Date(to[per]) )
selectU = selectU/window(S, start = as.Date(from[per]) , end = as.Date(to[per]) );
selectU = clean.boudt2(selectU , alpha = alpha )[[1]];
in.sample.T = nrow(selectU);
Rcor = matrix(rep(0,cAssets^2),nrow=cAssets,ncol=cAssets)
M3 = matrix(rep(0,cAssets^3),nrow=cAssets,ncol=cAssets^2)
M4 = matrix(rep(0,cAssets^4),nrow=cAssets,ncol=cAssets^3)
for(t in c(1:in.sample.T)){
centret = as.vector(selectU[t,]);
Rcor = Rcor + centret%*%t(centret)
M3 = M3 + ( centret%*%t(centret) )%x%t(centret)
M4 = M4 + ( centret%*%t(centret) )%x%t(centret)%x%t(centret)
}
resSigma = Rcor = 1/in.sample.T*Rcor
sdzero = function(data){ return( sqrt(mean(data^2)) ) }
D = diag( apply(selectU,2,'sdzero')^(-1),ncol=cAssets )
Rcor = D%*%Rcor%*%D
M3 = 1/in.sample.T*M3
M4 = 1/in.sample.T*M4
# we only need mean and conditional covariance matrix of last observation
mu = matrix(tail(M,n=1),ncol=1 ) ;
D = diag( as.vector(as.vector(tail(S,n=1) ) ),ncol=cAssets )
sigma = D%*%Rcor%*%D
switch( percriskcontribcriterion ,
StdDev = { percriskcontrib = function(w){ cont = Portsd(w,mu=mu,sigma=sigma)[[3]] ; return( cont ) }},
GVaR = {percriskcontrib = function(w){ cont = PortgausVaR(w,alpha=alpha,mu=mu,sigma=sigma)[[3]] ; return( cont ) }},
GES = {percriskcontrib = function(w){ cont = PortgausES(w,mu=mu,alpha=alpha,sigma=sigma)[[3]] ; return( cont ) }},
mVaR = {percriskcontrib = function(w){ cont = resPortMVaR(w,mu=mu,alpha=alpha,sigma=sigma,resSigma=resSigma,M3=M3,M4=M4)[[3]] ; return( cont ) }},
mES = {percriskcontrib = function(w){ cont = resoperPortMES(w,mu=mu,alpha=alpha,sigma=sigma,resSigma=resSigma,M3=M3,M4=M4)[[3]] ; return( cont ) }}
)
RCout[ per , ] = as.vector(percriskcontrib(wper))
} # end loop periods
#write.table( RCout , paste("riskcont/", percriskcontribcriterion,"/",datacase,"/",frequency,"/EW.csv" , sep=""),
# append = FALSE, quote = TRUE, sep = ",", eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE, qmethod = "escape")
write.table( RCout ,file = paste("C:\\Documents and Settings\\Administrator\\Desktop\\risk budget programs\\riskcont\\mES\\equitybondscommodity\\quarterly\\unconstrained\\", "mES.4yr-InfInf" , ".csv" , sep="") );
append = FALSE, quote = TRUE, sep = ",", eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE, qmethod = "escape")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.