User:Xyxht/sandbox


 * 1) set directory, import data
 * 2) please check data structure before run below codes when importing
 * 3) SPLus results
 * 4) There is a potential issue of PCA analysis. the length of the lookback period
 * 5) cannot be greater than the length of NotNAs stocks
 * 6)  Wei Lu
 * 7)  20/04/2015
 * 1) cannot be greater than the length of NotNAs stocks
 * 2)  Wei Lu
 * 3)  20/04/2015
 * 1)  Wei Lu
 * 2)  20/04/2015
 * 1)  Wei Lu
 * 2)  20/04/2015
 * 1)  Wei Lu
 * 2)  20/04/2015


 * 1) Lots value R resources available in this website. I've chosen small
 * 2) samples to test the relevant functions already before applying them
 * 3) to our main dataset. Also, I've gone through the details of each function
 * 4) So I can confirm these functions correct.
 * 5) Below step enables to connect function packages called "sit.gz"
 * 1) Below step enables to connect function packages called "sit.gz"
 * 1) Below step enables to connect function packages called "sit.gz"

con = gzcon(file('O:/Quant Research/R/Wei/external package/sit.gz', 'rb')) source(con) close(con) load.packages('ineq,quantmod,quadprog,corpcor,lpSolve,timeSeries,fPortfolio,Rdonlp2,tawny,PerformanceAnalytics,xlsx,clusterSim')
 * 1) https://github.com/systematicinvestor/SIT/raw/master/sit.gz - Systematic Investor Toolbox
 * 1) https://github.com/systematicinvestor/SIT/raw/master/sit.gz - Systematic Investor Toolbox


 * 1) con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
 * 2) source(con)
 * 3) close(con)


 * 1) Load all packages/libraries from R
 * 2) "load.packages" is an efficient function to download more than one packages
 * 3)  from R--from the above connection source Be careful : you must install all
 * 4)  packages first
 * 5)  click "Packages" in the right hand side, and then click "Install Packages".
 * 6)  The way that R works, need you install the relevant packages you want before running codes,
 * 7)  Normally I put all packages in the front, it's easy to check which one I missed.
 * 1)  The way that R works, need you install the relevant packages you want before running codes,
 * 2)  Normally I put all packages in the front, it's easy to check which one I missed.


 * 1)  The original R only includes basic packages, you can update packages anytime (someone may
 * 2)  add new packages in, you just share them as a free lunch).

load.packages('ineq,quantmod,quadprog,corpcor,lpSolve,timeSeries,fPortfolio,Rdonlp2,tawny,PerformanceAnalytics,xlsx,clusterSim')


 * 1) Set working directory - in order to get data
 * 1) Set working directory - in order to get data

setwd("O:/Quant Research/Strategic Work/Smart Beta/Data/Data") tcost = read.csv("Tcost Input Sheet.csv")
 * 1) Import transaction cost form
 * 1) tcost = read.csv("New Tcost Input Sheet.csv")

fildir = "O:/Quant Research/Strategic Work/Smart Beta/Data/Data/aisa index/Update 20140402/"
 * 1) Import rets series back with extra 5 years
 * 1) Import rets series back with extra 5 years
 * 1) fildir = "O:/Quant Research/Strategic Work/Smart Beta/Data/Data/aisa index/Weekly New/"

weightfile = "H:/Quant Research/Strategic Work/SmartBeta/final results/Weekly updated/PCA Portfolio weights/" PCAfile = "H:/Quant Research/Strategic Work/SmartBeta/final results/Weekly updated/PCA/" addMeanSemiVar = "H:/Quant Research/Strategic Work/SmartBeta/final results/Weekly updated/Universe stocks Summary/" portRet = "H:/Quant Research/Strategic Work/SmartBeta/final results/Weekly updated/PCA Ret/" stockName = "H:/Quant Research/Strategic Work/SmartBeta/final results/Weekly updated/Universe Stock Name/" ginicoef = "H:/Quant Research/Strategic Work/SmartBeta/final results/Weekly updated/PCA Concentration Ratio/"
 * 1) Export directory
 * 1) Export directory

index.name = c("AUSTRALIA","CHINA","Hong_Kong","INDIA","INDONESIA","MALAYSIA","PHILIPPINES",              "SINGAPORE","SOUTH_KOREA","TAIWAN","THAILAND") country.constraints = c(0.05,0.1,0.03,0.07,0.07,0.05,0.1,0.05,0.04,0.04,0.06)
 * 1) Load original dataset (price information) from splus output
 * 1) Load original dataset (price information) from splus output

setwd(fildir)
 * 1) Set working directoy
 * 1) Set working directoy

Myfiles <- list.files Myfiles = Myfiles[1:length(Myfiles)] # exclude new folder called "add5years"
 * 1) list all files in the working directory

a <- 1:length(Myfiles) b4 <- a[seq(1, length(a), 4)] b3 = c(b4-1,length(Myfiles)) b3 = b3[2:length(b3)]
 * 1) load data

optimize.portfolio.nlp = function (ia, constraints, fn, nl.constraints = NULL,                                   direction = 'min',full.solution = F) { if( direction == 'min' ) fnscale = 1 else fnscale = -1 if( as.numeric(sessionInfo$R.version$minor ) < 9 ) { cntl <- donlp2Control(silent = T, fnscale = fnscale, iterma =10000, nstep = 100, epsx = 1e-10) } else { cntl <- donlp2Control cntl$silent = T   cntl$fnscale = fnscale cntl$iterma =10000 cntl$nstep = 100 cntl$epsx = 1e-10 } par.l = constraints$lb par.u = constraints$ub p = rep(1, nrow(constraints$A)) if(!is.null(constraints$x0)) { if( sum(is.na(constraints$x0)) == 0) p = constraints$x0 } A = t(constraints$A) lin.l = constraints$b lin.u = constraints$b lin.u[ -c(1:constraints$meq) ] = +Inf x = NA if( !is.null(nl.constraints) ) { sol = donlp2(p, fn,                par.lower=par.l, par.upper=par.u,                 A=A, lin.u=lin.u, lin.l=lin.l,                 control=cntl,                 nlin=nl.constraints$constraints,                 nlin.upper=nl.constraints$upper, nlin.lower=nl.constraints$lower    ) } else { sol = donlp2(p, fn,                par.lower=par.l, par.upper=par.u,                 A=A, lin.u=lin.u, lin.l=lin.l,                 control=cntl) } if(!inherits(sol, 'try-error')) { x = sol$par } if( full.solution ) x = sol return( x ) }
 * 1) function of optimize.portfolio.nlp  -- constrained nonlinear - for equal risk contribution method
 * 2) this function should correct, I've tested using small samples.
 * 1) this function should correct, I've tested using small samples.



create.historical.ia = function ( hist.returns,  annual.factor,  symbols = colnames(hist.returns),  symbol.names = symbols) { ia = list ia$n = len(symbols) ia$annual.factor = annual.factor ia$symbols = symbols ia$symbol.names = symbol.names ia$hist.returns = hist.returns ia$arithmetic.return = apply(hist.returns, 2, mean, na.rm = T) #ia$geometric.return = apply(hist.returns, 2, function(x) prod(1+x)^(1/len(x))-1 ) ia$risk = apply(hist.returns, 2, sd, na.rm = T) #ia$correlation = cor(hist.returns, use = 'complete.obs', method = 'pearson') ia$arithmetic.return = (1 + ia$arithmetic.return)^ia$annual.factor - 1 #ia$geometric.return = (1 + ia$geometric.return)^ia$annual.factor - 1 #ia$risk = sqrt(ia$annual.factor) * ia$risk #ia$risk = iif(ia$risk == 0, 0.000001, ia$risk) #ia$cov = ia$correlation * (ia$risk %*% t(ia$risk)) ia$expected.return = ia$arithmetic.return return(ia) }

createPCfromRetProd = function(retMatrix, cutOff = 0.01) {  # Names of all stocks n = dimnames(retMatrix)2 # No of NA returns per stock nCountNA = apply(is.na(retMatrix), 2, sum) nOK = n[nCountNA == 0] nNA = n[nCountNA > 0] # Get priciple component model for stocks without any NA returns modPC = prcomp(retMatrix[, nOK, drop = F]) # Variance of each component mPCVar = modPC$sdev^2 # Components that explain at least 90% of total variance nPC = sum(mPCVar/sum(mPCVar) >= cutOff) # Factor loadings ldgPC = matrix(NA, length(n), nPC) dimnames(ldgPC) = list(n, paste("comp", 1:nPC, sep=".")) #ldgPC[nOK,] = modPC$loadings[, seq(nPC), drop = F]  ldgPC[nOK,] = modPC$rotation[, seq(nPC), drop = F]     # Factor Returns retPC = retMatrix[, nOK, drop = F] %*% ldgPC[nOK,] # Identify NA stocks with limited history, but still estimateable factor loadings nCount = apply(!is.na(retMatrix[, nNA]), 2, sum) # Factor Covariance varPC = var(retPC) # Stock Specific Variance varResid = colVars(retMatrix - retPC %*% t(ldgPC), na.rm = T)  if(sum(is.na(varResid)) > 0) varResid[is.na(varResid)] = median(varResid, na.rm = T)  list(loadings = ldgPC, covarPC = varPC, varResid = varResid, returnPC = retPC) }

all.index.price = list temp.rets = list weekly.country.weight = list weekly.country.invert.weight = list
 * 1) Calculate stocks weights of each country based on different index strategy, the order is
 * 2) for each market, apply all index methodologies, once done, then go to next market.
 * 1) for each market, apply all index methodologies, once done, then go to next market.
 * 1) Chose market
 * 1) Chose market

for(k in 1:length(index.name)){ #k = 1 SharpeRatio.weights = list SharpeRatio.1stHlf = list SharpeRatio.2ndHlf = list #for(k in 1:length(index.name)){ ind.index = k constrain = country.constraints[k] # Prepare Data ## port return temp.retsind.index = read.csv(paste(fildir,"/",index.namek,"_Data_w tr.csv",sep = "")) ind.splus = b4ind.index:b3ind.index all.index.priceind.index = lapply(1:length(ind.splus),function (x) read.csv(Myfilesind.splus[[x]])) d = all.index.priceind.index1[,1] # on each week, replace universe.dummy with 0 to NA in market value and weekly return for(j in 1:length(d)){ uni.dummy = which(all.index.priceind.index4[j,] == 0) all.index.priceind.index1[j,uni.dummy] = NA # if universe.dummy is o, then setting size = NA } for(j in 1:length(d)){ uni.dummy = which(all.index.priceind.index4[j,] == 0) all.index.priceind.index3[j,uni.dummy] = NA # if universe.dummy is o, then setting size = NA } #  dd = which(temp.retsind.index[,1]==d1) #  temp.retsind.index[(dd-52*5):dd,which(all.index.priceind.index4[1,] == 0)] = NA  # get size and rets ready to exlude stocks not in the index - weekly basis size = all.index.priceind.index1 ret = all.index.priceind.index3 ############################ Capped for Hk total return ############################ if(ind.index == 3){ #ret1 = ret for(t in 1:nrow(ret)){ tret = ret[t,2:ncol(ret)] nna.ind = which(!is.na(tret)==TRUE) ttret = as.matrix(tret[nna.ind]) ttret = pmax(ttret, -0.5) ttret = pmin(ttret, 1) ret[t,nna.ind+1] = as.numeric(as.matrix(ttret)) }   #ret1[,1] = paste(substr(ret1[,1],1,4),"-",substr(ret1[,1],5,6),"-",substr(ret1[,1],7,8),sep = "") n = ncol(ret) } ###################################################################  #check how many stock in the universe #sapply(1:len(ret[,1]),function(x){len(which(!is.na(ret[x,])))}) # convert size and ret to time series structure size[,1] = paste(substr(size[,1],1,4),"-",substr(size[,1],5,6),"-",substr(size[,1],7,8),sep = "") n = ncol(ret) sizes = timeSeries(size[,2:length(size[1,])],size[,1]) ret[,1] = paste(substr(ret[,1],1,4),"-",substr(ret[,1],5,6),"-",substr(ret[,1],7,8),sep = "") n = ncol(ret) rets = timeSeries(ret[,2:length(ret[1,])],ret[,1]) rets = rets[1:length(rets[,1]),] #==================================================================================================================== # Create Portfolios #==================================================================================================================== lookback.period = 52*2 # get last five years data to calculate covariance matrix, 52 weeks, 5 years #ret = prices / mlag(prices) - 1 # weely return start.i = lookback.period + 1 # Rebalance quarterly rebalance.horizon = 52/4 all.rebal.period = (start.i):length(rets[,1]) effect.rebal.period = all.rebal.period[seq(1, length(all.rebal.period), rebalance.horizon)] ## Create weight variable weight = NA * rets weights = list weights$equal.weight = weight    # equal weighted weights$min.var = weight         # minimum variance weights$max.div = weight         # maximum diversification weights$erc = weight             # equal risk contribution weights$min.downside = weight    # minimum semi-variance weights$inv.volatility = weight  # Inverse volatility weights$min.corr = weight        # minimum correlation weights$min.corr2 = weight       # minimum correlation2 weights$cap = weight             # cap weighted for( j in 1:length(effect.rebal.period) ) { # calculate portfolio weights on each index method m = effect.rebal.periodj #This is based on weekly returns - lookback period to get covariance matrix - 52*5 hist.temp = rets[(m-lookback.period+1):m, ] size.temp = sizes[(m-lookback.period+1):m, ] ### exclude stocks with all NA returns in a period notna.index = which(!is.na(apply(coredata(hist.temp),2,sum))) hhist = hist.temp[,notna.index] hist = hhist hist = hist[,which(colSums(hist)!=0)] notna.index = notna.index[which(colSums(hist)!=0)] ###exclude stocsk with all NA sizes in a period notna.size.index = which(!is.na(apply(coredata(size.temp),2,sum))) ssize = size.temp[,notna.size.index] ssizes = ssize ssizes = ssizes[,which(colSums(ssizes)!=0)] notna.size.index = notna.size.index[which(colSums(ssizes)!=0)] notna.index = intersect(notna.index,notna.size.index) hist = hist.temp[,notna.index] # PCA Create the Variance-Covariance matrix pc = createPCfromRetProd(hist) cov.temp = pc$loadings %*% pc$covarPC %*% t(pc$loadings) + diag(pc$varResid) ia = create.historical.ia(hist,52) ia$cov = cov.temp ia$correlation = cov2cor(ia$cov) #s0 = apply(coredata(hist),2,sd) #  covar =  cov_shrink(as.zoo(hist)) #  ia$cov = covar#1 # set constrains 1. no short selling 2. sum = 1 # 0 <= x.i <= 1 constraints = new.constraints(ia$n, lb = 0, ub = constrain) constraints = add.constraints(diag(ia$n), type='>=', b = 0, constraints) constraints = add.constraints(diag(ia$n), type='<=', b = 1, constraints) # SUM x.i = 1 constraints = add.constraints(rep(1, ia$n), type = '=', b = 1, constraints) # set nonlinear constraints # Equal weighetd weights$equal.weight[m,notna.index] = 1/len(notna.index) # Minimum variacne - two ways to get it     # 1. Based on bulit in function from website, called min.risk.portfolio weights$min.var[m,notna.index] = min.risk.portfolio(ia, constraints) # Minimum correlation weights$min.corr[m,notna.index] = min.corr.portfolio(ia, constraints) # Minimum correlation weights$min.corr2[m,notna.index] = min.corr2.portfolio(ia, constraints) # 2. use solve.QP to get minimum variance #     N = ncol(hist) #     zeros = array(0, dim = c(N,1)) #     aMat = cbind(t(array(1, dim = c(1,N))), diag(N)) #short selling constraint (i.e. non-negative weights): #     b0 = as.matrix(c(1, rep.int(0,N))) #     res = solve.QP(ia$cov, zeros, aMat, bvec=b0, meq = 1) #     weights$min.var[m,notna.index] = res$solution # Maximum diversified - "max.div.portfolio" function is built in, it is based on "min.var.portfolio" #    max.div.portfolio = function #    (    #       ia,    #       constraints    #     ) #    {    #       risk.index = get.risky.asset.index(ia) #      x = min.var.portfolio(ia, constraints, ia$cov) #      x = x[risk.index] / ia$risk[risk.index] #      set.risky.asset(x / sum(x), risk.index) #    }    weights$max.div[m,notna.index] = max.div.portfolio(ia,constraints) # Equal risk contribution - this function "optimize.portfolio.nlp" is bulit in    weights$erc[m,notna.index] = find.erc.portfolio(ia, constraints) # Mean-Semi Variance portfolio weights$min.downside[m,notna.index] = min.risk.downside.portfolio(ia,constraints) # Inverse Volatility volatility = ia$hist volatility[] = NA    d = dimnames(hist)1 d2 = d    for(i in d[-(1:(52*1))]) { d1 = d2[which(d2 <= i)] d1 = rev(rev(d1)[1:(52*1)]) volatility[which(d==i),] = colStdevs(ia$hist[which(d%in%d1), ], na.rm=T) NULL }    volatility = volatility * sqrt(52*1) volatility[1:(52*1),] = t(matrix(volatility[52*1+1,], dim(ia$hist)[2], 52*1)) volatility[!is.finite(volatility)] = NA   volatility[is.finite(volatility) & volatility == 0] = NA    new.na = which(!is.na(volatility[nrow(volatility),])) new.notna.index = notna.index[new.na] weights$inv.volatility[m,new.notna.index] = 1/volatility[nrow(volatility),!is.na(volatility[nrow(volatility),])] weights$inv.volatility[m,new.notna.index] = weights$inv.volatility[m,new.notna.index]/sum(weights$inv.volatility[m,new.notna.index]) # Capital weighted - based on market value from Factset for universe dummy for(s in 1:length(notna.index)){ weights$cap[m,notna.indexs] = sizes[m,notna.indexs]/sum(sizes[m,notna.index]) }      } # all different weights done !!! #======================================================================== # Adjust for Max div #======================================================================== single.weights=list single.weights1 = weights$max.div names(single.weights)= index.name[k] cap.max.div = list for(m in 1:1){ max.cap = country.constraintsm temp.max.div = single.weightsm temp.max.div[temp.max.div<=0.0001] = 0 temp.max.div[is.na(temp.max.div)] = 0 new.max.div = temp.max.div new.max.div[]=0 max.weight = apply(temp.max.div,1,max) for(i in 1:length(max.weight)){ ind.non.zero = which(temp.max.div[i,]!=0) order.non.zero = order(temp.max.div[i,],decreasing = TRUE) if(max(temp.max.div[i,])>max.cap){ if((1/length(ind.non.zero)) < max.cap){ #extra.cap[] = NA          extra.cap = temp.max.div[i,order.non.zero] - max.cap p.extra.cap = which(extra.cap>= 0) n.extra.cap = which(extra.cap<0 & extra.cap >-max.cap) #tp[] = NA          tp = temp.max.div[i,order.non.zero] tp[p.extra.cap] = max.cap tp[n.extra.cap] = tp[n.extra.cap] + sum(extra.cap[p.extra.cap])/length(n.extra.cap) name.tp = dimnames(tp)2 name.old = dimnames(temp.max.div)2 new.max.div[i,] = as.numeric(as.matrix(tp[match(name.old,name.tp)])) }else if((1/length(ind.non.zero)) >= max.cap){ # if 1/n is greater than max.cap,then do equal weighting new.max.div[i,ind.non.zero] = 1/length(ind.non.zero) }     }else { new.max.div[i,] = temp.max.div[i,] }   }    cap.max.divm = new.max.div } single.weights = cap.max.div names(single.weights)= index.namek new.cap.max.div = list for(m in 1:1){ max.cap = country.constraintsm temp.max.div = single.weightsm temp.max.div[temp.max.div<=0.0001] = 0 temp.max.div[is.na(temp.max.div)] = 0 new.max.div = temp.max.div new.max.div[]=0 max.weight = apply(temp.max.div,1,max) for(i in 1:length(max.weight)){ ind.non.zero = which(temp.max.div[i,]!=0) order.non.zero = order(temp.max.div[i,],decreasing = TRUE) if(max(temp.max.div[i,])>max.cap){ if((1/length(ind.non.zero)) < max.cap){ #extra.cap[] = NA          extra.cap = temp.max.div[i,order.non.zero] - max.cap p.extra.cap = which(extra.cap>= 0) n.extra.cap = which(extra.cap<0 & extra.cap >-max.cap) #tp[] = NA          tp = temp.max.div[i,order.non.zero] tp[p.extra.cap] = max.cap tp[n.extra.cap] = tp[n.extra.cap] + sum(extra.cap[p.extra.cap])/length(n.extra.cap) name.tp = dimnames(tp)2 name.old = dimnames(temp.max.div)2 new.max.div[i,] = as.numeric(as.matrix(tp[match(name.old,name.tp)])) }else if((1/length(ind.non.zero)) >= max.cap){ # if 1/n is greater than max.cap,then do equal weighting new.max.div[i,ind.non.zero] = 1/length(ind.non.zero) }     }else { new.max.div[i,] = temp.max.div[i,] }   }    new.cap.max.divm = new.max.div } single.weights = new.cap.max.div names(single.weights)= index.namek new1.cap.max.div = list for(m in 1:1){ max.cap = country.constraintsm temp.max.div = single.weightsm temp.max.div[temp.max.div<=0.0001] = 0 temp.max.div[is.na(temp.max.div)] = 0 new.max.div = temp.max.div new.max.div[]=0 max.weight = apply(temp.max.div,1,max) for(i in 1:length(max.weight)){ ind.non.zero = which(temp.max.div[i,]!=0) order.non.zero = order(temp.max.div[i,],decreasing = TRUE) if(max(temp.max.div[i,])>max.cap){ if((1/length(ind.non.zero)) < max.cap){ #extra.cap[] = NA          extra.cap = temp.max.div[i,order.non.zero] - max.cap p.extra.cap = which(extra.cap>= 0) n.extra.cap = which(extra.cap<0 & extra.cap >-max.cap) #tp[] = NA          tp = temp.max.div[i,order.non.zero] tp[p.extra.cap] = max.cap tp[n.extra.cap] = tp[n.extra.cap] + sum(extra.cap[p.extra.cap])/length(n.extra.cap) name.tp = dimnames(tp)2 name.old = dimnames(temp.max.div)2 new.max.div[i,] = as.numeric(as.matrix(tp[match(name.old,name.tp)])) }else if((1/length(ind.non.zero)) >= max.cap){ # if 1/n is greater than max.cap,then do equal weighting new.max.div[i,ind.non.zero] = 1/length(ind.non.zero) }     }else { new.max.div[i,] = temp.max.div[i,] }   }    new1.cap.max.divm = new.max.div } single.weights = new1.cap.max.div names(single.weights)= index.namek new2.cap.max.div = list for(m in 1:1){ max.cap = country.constraintsm temp.max.div = single.weightsm temp.max.div[temp.max.div<=0.0001] = 0 temp.max.div[is.na(temp.max.div)] = 0 new.max.div = temp.max.div new.max.div[]=0 max.weight = apply(temp.max.div,1,max) for(i in 1:length(max.weight)){ ind.non.zero = which(temp.max.div[i,]!=0) order.non.zero = order(temp.max.div[i,],decreasing = TRUE) if(max(temp.max.div[i,])>max.cap){ if((1/length(ind.non.zero)) < max.cap){ #extra.cap[] = NA          extra.cap = temp.max.div[i,order.non.zero] - max.cap p.extra.cap = which(extra.cap>= 0) n.extra.cap = which(extra.cap<0 & extra.cap >-max.cap) #tp[] = NA          tp = temp.max.div[i,order.non.zero] tp[p.extra.cap] = max.cap tp[n.extra.cap] = tp[n.extra.cap] + sum(extra.cap[p.extra.cap])/length(n.extra.cap) name.tp = dimnames(tp)2 name.old = dimnames(temp.max.div)2 new.max.div[i,] = as.numeric(as.matrix(tp[match(name.old,name.tp)])) }else if((1/length(ind.non.zero)) >= max.cap){ # if 1/n is greater than max.cap,then do equal weighting new.max.div[i,ind.non.zero] = 1/length(ind.non.zero) }     }else { new.max.div[i,] = temp.max.div[i,] }   }    new2.cap.max.divm = new.max.div } single.weights = new2.cap.max.div names(single.weights)= index.namek new3.cap.max.div = list for(m in 1:1){ max.cap = country.constraintsm temp.max.div = single.weightsm temp.max.div[temp.max.div<=0.0001] = 0 temp.max.div[is.na(temp.max.div)] = 0 new.max.div = temp.max.div new.max.div[]=0 max.weight = apply(temp.max.div,1,max) for(i in 1:length(max.weight)){ ind.non.zero = which(temp.max.div[i,]!=0) order.non.zero = order(temp.max.div[i,],decreasing = TRUE) if(max(temp.max.div[i,])>max.cap){ if((1/length(ind.non.zero)) < max.cap){ #extra.cap[] = NA          extra.cap = temp.max.div[i,order.non.zero] - max.cap p.extra.cap = which(extra.cap>= 0) n.extra.cap = which(extra.cap<0 & extra.cap >-max.cap) #tp[] = NA          tp = temp.max.div[i,order.non.zero] tp[p.extra.cap] = max.cap tp[n.extra.cap] = tp[n.extra.cap] + sum(extra.cap[p.extra.cap])/length(n.extra.cap) name.tp = dimnames(tp)2 name.old = dimnames(temp.max.div)2 new.max.div[i,] = as.numeric(as.matrix(tp[match(name.old,name.tp)])) }else if((1/length(ind.non.zero)) >= max.cap){ # if 1/n is greater than max.cap,then do equal weighting new.max.div[i,ind.non.zero] = 1/length(ind.non.zero) }     }else { new.max.div[i,] = temp.max.div[i,] }   }    new3.cap.max.divm = new.max.div } single.weights = new3.cap.max.div names(single.weights)= index.namek new4.cap.max.div = list for(m in 1:1){ max.cap = country.constraintsm temp.max.div = single.weightsm temp.max.div[temp.max.div<=0.0001] = 0 temp.max.div[is.na(temp.max.div)] = 0 new.max.div = temp.max.div new.max.div[]=0 max.weight = apply(temp.max.div,1,max) for(i in 1:length(max.weight)){ ind.non.zero = which(temp.max.div[i,]!=0) order.non.zero = order(temp.max.div[i,],decreasing = TRUE) if(max(temp.max.div[i,])>max.cap){ if((1/length(ind.non.zero)) < max.cap){ #extra.cap[] = NA          extra.cap = temp.max.div[i,order.non.zero] - max.cap p.extra.cap = which(extra.cap>= 0) n.extra.cap = which(extra.cap<0 & extra.cap >-max.cap) #tp[] = NA          tp = temp.max.div[i,order.non.zero] tp[p.extra.cap] = max.cap tp[n.extra.cap] = tp[n.extra.cap] + sum(extra.cap[p.extra.cap])/length(n.extra.cap) name.tp = dimnames(tp)2 name.old = dimnames(temp.max.div)2 new.max.div[i,] = as.numeric(as.matrix(tp[match(name.old,name.tp)])) }else if((1/length(ind.non.zero)) >= max.cap){ # if 1/n is greater than max.cap,then do equal weighting new.max.div[i,ind.non.zero] = 1/length(ind.non.zero) }     }else { new.max.div[i,] = temp.max.div[i,] }   }    new4.cap.max.divm = new.max.div } single.weights = new4.cap.max.div names(single.weights)= index.namek new5.cap.max.div = list for(m in 1:1){ max.cap = country.constraintsm temp.max.div = single.weightsm temp.max.div[temp.max.div<=0.0001] = 0 temp.max.div[is.na(temp.max.div)] = 0 new.max.div = temp.max.div new.max.div[]=0 max.weight = apply(temp.max.div,1,max) for(i in 1:length(max.weight)){ ind.non.zero = which(temp.max.div[i,]!=0) order.non.zero = order(temp.max.div[i,],decreasing = TRUE) if(max(temp.max.div[i,])>max.cap){ if((1/length(ind.non.zero)) < max.cap){ #extra.cap[] = NA          extra.cap = temp.max.div[i,order.non.zero] - max.cap p.extra.cap = which(extra.cap>= 0) n.extra.cap = which(extra.cap<0 & extra.cap >-max.cap) #tp[] = NA          tp = temp.max.div[i,order.non.zero] tp[p.extra.cap] = max.cap tp[n.extra.cap] = tp[n.extra.cap] + sum(extra.cap[p.extra.cap])/length(n.extra.cap) name.tp = dimnames(tp)2 name.old = dimnames(temp.max.div)2 new.max.div[i,] = as.numeric(as.matrix(tp[match(name.old,name.tp)])) }else if((1/length(ind.non.zero)) >= max.cap){ # if 1/n is greater than max.cap,then do equal weighting new.max.div[i,ind.non.zero] = 1/length(ind.non.zero) }     }else { new.max.div[i,] = temp.max.div[i,] }   }    new5.cap.max.divm = new.max.div } single.weights = new5.cap.max.div names(single.weights)= index.namek new6.cap.max.div = list for(m in 1:1){ max.cap = country.constraintsm temp.max.div = single.weightsm temp.max.div[temp.max.div<=0.0001] = 0 temp.max.div[is.na(temp.max.div)] = 0 new.max.div = temp.max.div new.max.div[]= 0 max.weight = apply(temp.max.div,1,max) for(i in 1:length(max.weight)){ ind.non.zero = which(temp.max.div[i,]!=0) order.non.zero = order(temp.max.div[i,],decreasing = TRUE) if(max(temp.max.div[i,])>max.cap){ if((1/length(ind.non.zero)) < max.cap){ #extra.cap[] = NA          extra.cap = temp.max.div[i,order.non.zero] - max.cap p.extra.cap = which(extra.cap>= 0) n.extra.cap = which(extra.cap<0 & extra.cap >-max.cap) #tp[] = NA          tp = temp.max.div[i,order.non.zero] tp[p.extra.cap] = max.cap tp[n.extra.cap] = tp[n.extra.cap] + sum(extra.cap[p.extra.cap])/length(n.extra.cap) name.tp = dimnames(tp)2 name.old = dimnames(temp.max.div)2 new.max.div[i,] = as.numeric(as.matrix(tp[match(name.old,name.tp)])) }else if((1/length(ind.non.zero)) >= max.cap){ # if 1/n is greater than max.cap,then do equal weighting new.max.div[i,ind.non.zero] = 1/length(ind.non.zero) }     }else { new.max.div[i,] = temp.max.div[i,] }   }    new6.cap.max.divm = new.max.div } #=============================================================================  temp.max.div = new6.cap.max.div1 #============================================================================= # Exclude weights of non-rebalance period #=============================================================================   act.weights = lapply(1:length(weights),function (x)     if(x!=3){      ty = weightsx[effect.rebal.period,]    }else{      ty = temp.max.div[effect.rebal.period,]    }) act.weights3[is.na(act.weights1)] = NA names(act.weights) = c("equal weighted", "min variance","max diversification",                         "equal risk","Mean-Lower Semi-Variance","Inverse Volatility",                         "min correlation","min correlation2","cap weighted") #============================================================================= # Compute inverse weights and max(r) - x  #============================================================================= diff.max.weights = lapply(2:length(weights),function (x){   sapply(1:nrow(act.weightsx),function(y){ ss = max(act.weightsx[y,],na.rm = T) - as.numeric(as.matrix(act.weightsx[y,])) ssum = sum(ss,na.rm = T)     diff.max = ss/ssum #diff.max = data.frame(diff.max) return(diff.max) }) })  diff.max.weights = lapply(1:length(diff.max.weights),function(x){    diff = t(diff.max.weightsx)    dimnames(diff)1 = rownames(act.weights1)    dimnames(diff)2 = colnames(act.weights1)    return(diff)  }) #==================================================================================================== # Export final.summary table, if you have not changed the above code, please comment on this #==================================================================================================== weights.name = c("EqualW","MinVar","MaxDiv","EqualR", "MeanLowerSemiVar","InverseVolatility",                   "MinCorr","MinCorr2", "CapW") sapply(1:len(weights),function(x){   write.xlsx(act.weightsx,paste(weightfile,                                      paste(index.nameind.index,weights.name[x],"port weights.xlsx"),sep = ""), row.names = T)}) diff.max.weights.name = c("Invert_MinVar","Invert_MaxDiv","Invert_EqualR",                            "Invert_MeanLowerSemiVar","Invert_InverseVolatility",                            "Invert_MinCorr","Invert_MinCorr2", "Invert_CapW") names(diff.max.weights) = diff.max.weights.name sapply(1:len(diff.max.weights),function(x){   write.xlsx(diff.max.weightsx,paste(weightfile,                                           paste(index.nameind.index, diff.max.weights.name[x],"weights.xlsx"),sep = ""), row.names = T)}) ############################# STORE WEEKLY WEIGHTS IN EACH COUNTRY  ################################### weekly.country.weightk = act.weights weekly.country.invert.weightk = diff.max.weights ######################################################################################################## #=================================================================  # Export Concentration Ratio - sum of squares and Gini Coefficient #================================================================= Gini1 = sapply(1:len(weights), function(x) sapply(1:len(act.weightsx[,1]), function(y) sum(act.weightsx[y,]^2,na.rm = T))) Gini2 = sapply(1:len(diff.max.weights), function(x) sapply(1:len(diff.max.weightsx[,1]), function(y) sum(diff.max.weightsx[y,]^2,na.rm = T))) dimnames(Gini1)1 = dimnames(act.weights1)1 dimnames(Gini1)2 = weights.name dimnames(Gini2)1 = dimnames(act.weights1)1 dimnames(Gini2)2 = diff.max.weights.name Gini = cbind(Gini1, Gini2) #export write.xlsx(Gini, paste(ginicoef,index.nameind.index," PCA Concentration Ratio.xlsx",sep = ""),row.names = T) #==================================================================================================== # Backtest - performance analysis #==================================================================================================== all.weights = list all.weights[1:length(act.weights)] = act.weights all.weights[(length(act.weights)+1):(length(act.weights)+length(diff.max.weights))] = diff.max.weights # Exclude NAs in weights matrix act.ind = which(rowSums(all.weights1,na.rm = T) == 1) # set initial weights factor actweightsFactor = list actweightsFactor = lapply(1:length(all.weights),function(x){   all.weightsx[act.ind1:length(all.weightsx[,1]),]}) names(all.weights) =c(names(act.weights),names(diff.max.weights)) # backtest function - modified based on Eben Van Wyk's version new.backtest.sub = function(p, lag = 0,rbalance = effect.rebal.period,                              rets, cost) { n = dimnames(rets)2 long.turnover = p[,1] long.turnover[] = NA   stock.rets = p    stock.rets[] = NA    for( i in 2:length(effect.rebal.period)){ # Turnover calculation of each weighting method new.p = as.matrix(unlist(sapply(1:length(n),function (x){       p[i-1,n[x]]*prod((1+rets[(effect.rebal.periodi-1+1):(effect.rebal.periodi-1),n[x]]))      }))) long.turnover[i]= colSums(t(abs(p[i-1,]- t(new.p))),na.rm = T)     # Get ready for stock rets in the performance calculation for (x in 1:length(n)){ stock.rets[i,x]= prod((1+rets[(effect.rebal.periodi-1+1):(effect.rebal.periodi),n[x]]))-1 }     }    #Adjust portfolio dates so that returns are correctly stated by realised date dimnames(p)1 = dimnames(stock.rets)1 # Long-only quarterly performance longP = p   longP[longP < 0] = 0 act.stock.rets = stock.rets[2:length(longP[,1]),] dimnames(act.stock.rets)1 = dimnames(longP)1[1:(length(longP[,1])-1)] long.rets = colSums(t(longP[1:(length(longP[,1])-1),] * act.stock.rets), na.rm = T)   long.turnover = long.turnover #Summary rowNames = c("long.rets") colNames = c("Return.pa", "Risk.pa", "Turnover.pa", "Cost.pa") summary = matrix(NA, length(rowNames), length(colNames), dimnames = list(rowNames, colNames)) # Return and risk for(i in rowNames) { summary[i, 1] = eval(parse(text = paste("mean(", i, ", na.rm = T) * 52/rebalance.horizon", sep = ""))) summary[i, 2] = eval(parse(text = paste("stdev(", i, ", na.rm = T) * sqrt(52/rebalance.horizon)", sep = ""))) }     # Turnover summary[i, 3] = mean(long.turnover, na.rm = T) * 52/rebalance.horizon # Cost summary[, 4] = summary[, 3] * cost list(summary = summary, long.rets = long.rets,        long.turnover = long.turnover,date = dimnames(p)1[2:length(p[,1])]) } # backtest results ### whole period # ret1 = all.index.priceind.index3 # #Cap weekly returns to +100% and -50% # for(t in 1:nrow(ret1)){ #  tret = ret1[t,2:ncol(ret1)] #  nna.ind = which(!is.na(tret)==TRUE) #  ttret = as.matrix(tret[nna.ind]) #  ttret = pmax(ttret, -0.5) #  ttret = pmin(ttret, 1) #  ret1[t,nna.ind+1] = as.numeric(as.matrix(ttret)) # } # ret1[,1] = paste(substr(ret1[,1],1,4),"-",substr(ret1[,1],5,6),"-",substr(ret1[,1],7,8),sep = "") # n = ncol(ret1) # rets = timeSeries(ret1[,2:length(ret1[1,])],ret1[,1]) # rets = rets[1:length(rets[,1]),] backtest.lag0 = lapply(1:length(all.weights),function(x){   new.backtest.sub(p = all.weightsx, rets = rets, cost = tcost[ind.index,3])  }) #================================================================================================================= # Get data ready for beta, sharpe ratio calculation #================================================================================================================= dates = backtest.lag01$date final.weights.rets = list final.weights.rets = lapply(1:length(all.weights),function(x){   ts(cumsum(backtest.lag0x$long.rets[(1:length(dates))]), start = c(as.numeric(substring(dates[1],1,4)),1),frequency = 52/rebalance.horizon )}) colNames = c("Beta","Sharp Ratio") summary.temp = matrix(0, length(colNames),length(names(all.weights))) dimnames(summary.temp) = list(colNames,names(all.weights)) weights.rets = cbind(backtest.lag01$long.rets,                      backtest.lag02$long.rets,                       backtest.lag03$long.rets,                       backtest.lag04$long.rets,                       backtest.lag05$long.rets,                       backtest.lag06$long.rets,                       backtest.lag07$long.rets,                       backtest.lag08$long.rets,                       backtest.lag09$long.rets,                       backtest.lag010$long.rets,                       backtest.lag011$long.rets,                       backtest.lag012$long.rets,                       backtest.lag013$long.rets,                       backtest.lag014$long.rets,                       backtest.lag015$long.rets,                       backtest.lag016$long.rets,                       backtest.lag017$long.rets  ) dimnames(weights.rets)1 = backtest.lag01$date dimnames(weights.rets)2 = names(all.weights) #======================================================================== # Beta - built in function in R - alternatively, get beta in Splus #======================================================================== beta.weights = list beta.weights = sapply(1:length(weights.rets[1,]),function(x){   CAPM.beta(weights.rets[,x],weights.rets[,9],Rf = 0)}) #======================================================================== # SharpeRatio - built in function R - have checked specifications in using this one #======================================================================== SharpeRatio.weightsk = sapply(1:(length(weights.rets[1,])),function(x){    SharpeRatio.annualized(weights.rets[,x],Rf = 0, scale = 52/rebalance.horizon,geometric = F)}) summary.temp[1,] = beta.weights summary.temp[2,] = SharpeRatio.weightsk # Construct final table including everything. summary.risk.weights = sapply(1:length(all.weights),function(x){   backtest.lag0x1[1,]}) final.summary = rbind(summary.risk.weights,summary.temp) colnames(final.summary) = names(all.weights) print(final.summary) # 1st and 2nd half Sharpe Ratio d.half = round(length(weights.rets[,1])/2) SharpeRatio.1stHlfk = sapply(1:(length(weights.rets[1,])),function(x){   SharpeRatio.annualized(weights.rets[1:d.half,x],Rf = 0, scale = 52/rebalance.horizon,geometric = F)}) SharpeRatio.2ndHlfk = sapply(1:(length(weights.rets[1,])),function(x){   SharpeRatio.annualized(weights.rets[(d.half+1):NROW(weights.rets),x],Rf = 0, scale = 52/rebalance.horizon,geometric = F)}) t.final.summary = rbind(final.summary,SharpeRatio.1stHlfk,SharpeRatio.2ndHlfk) rownames(t.final.summary)[7] = c("Sharpe 1st half") rownames(t.final.summary)[8] = c("Sharpe 2nd half") #=========================================================================================================== #Export results to PCAfiles #=========================================================================================================== write.xlsx(t.final.summary,file = paste(addMeanSemiVar, paste(index.nameind.index,"PCA all weights summary statistics.xlsx"), sep = ""),row.names = T) final.cum.ret.weights = do.call(cbind,final.weights.rets) final.cum.ret.weights <- cbind(   as.matrix(final.weights.rets1),as.matrix(final.weights.rets2),    as.matrix(final.weights.rets3),as.matrix(final.weights.rets4),    as.matrix(final.weights.rets5),as.matrix(final.weights.rets6),    as.matrix(final.weights.rets7),as.matrix(final.weights.rets8),    as.matrix(final.weights.rets9),as.matrix(final.weights.rets10),    as.matrix(final.weights.rets11),as.matrix(final.weights.rets12),    as.matrix(final.weights.rets13),as.matrix(final.weights.rets14),    as.matrix(final.weights.rets15),as.matrix(final.weights.rets16),    as.matrix(final.weights.rets17)  ) dimnames(final.cum.ret.weights )1 = backtest.lag01$date dimnames(final.cum.ret.weights )2 = names(all.weights) # Export final.cum.ret, if you have not changed the above code, please comment on this write.xlsx(final.cum.ret.weights,paste(addMeanSemiVar,paste(index.nameind.index,                                                              "PCA Cum Ret.xlsx"),sep = ""),row.names = T)  write.xlsx(weights.rets,paste(portRet,paste("format",                                              index.nameind.index, "PCA Port Ret.xlsx"),sep = ""), row.names = T) }
 * 1) Export Stock Names in each market
 * 2) sapply(1:len(weights),function(x){
 * 3)   Name= dimnames(act.weightsx)2
 * 4)   write.xlsx(t(Name),paste(stockName,paste(
 * 5)     index.nameind.index, "Stock Name.xlsx"),sep = ""),row.names = F)})
 * 1)   write.xlsx(t(Name),paste(stockName,paste(
 * 2)     index.nameind.index, "Stock Name.xlsx"),sep = ""),row.names = F)})