# Commands that generated the results of Section 2.3 of
# "The Portfolio Probe User's Manual"
#
# These commands demonstrate some ways of testing the quality of risk models
# using random portfolios
#
# This depends on functions listed in 'pprobe_functions01.R'
# The data are created in 'pprobe_R_data01.txt'
#
# The functions that create the graphs are in 'pprobe_graphFun01.R'
# If you want hardcopies of the graphs, then you need to start a
# graphics device (like 'pdf') and do 'dev.off()' after the graphs
# are created.




# things to CHANGE, possibly

PP_LOC <- ".." 





# initial setup

if(!exists("us.pricemat")) attach(paste(PP_LOC, "us.pricemat.rda", sep="/"))
if(!exists("us.macdmat")) attach(paste(PP_LOC, "us.macdmat.rda", sep="/"))

quarter.ends <- c(251, 312, 375, 438, 502, 563, 627, 691, 755)
names(quarter.ends) <- paste(rep(2006:2008, c(1,4,4)), 'Q', c(4,1:3), sep="")

empty8 <- vector("list", 8)
names(empty8) <- tail(names(quarter.ends), 8)

# create return matrix
us.retmat <- diff(log(us.pricemat))

# create variance matrices

require(BurStFin)

varshr.Q <- empty8
varfac.Q <- empty8

for(i in 1:8) {
	varshr.Q[[i]] <- var.shrink.eqcor(us.retmat[seq(to=quarter.ends[i]-1, 
		length=250),])
	varfac.Q[[i]] <- factor.model.stat(us.retmat[seq(to=quarter.ends[i]-1, 
		length=250),])
}

# create lists of random portfolios at each quarter
rp30s15.Q <- empty8
for(i in 1:8) {
	rp30s15.Q[[i]] <- random.portfolio(1000, us.pricemat[quarter.ends[i], ], 
		gross.value=1e6, long.only=TRUE, port.size=30, 
		max.weight=.05, sum.weight=c("4"=.15))
}

# find ex ante volatilities of random portfolios

rp30s15.facexante.Q <- empty8
rp30s15.shrexante.Q <- empty8

for(i in 1:8) {

	# there is a subtle but important issue here
	# the prices should be given explicitly because the prices were
	# originally: us.pricemat[quarter.ends[i], ]
	# we don't want to depend on the loop 'i' here being the same
	# as the original 'i'

	rp30s15.facexante.Q[[i]] <- unlist(randport.eval(rp30s15.Q[[i]], 
		additional.args=list(prices=us.pricemat[quarter.ends[i], ], 
		variance=varfac.Q[[i]]), keep="var.values"), use.names=FALSE)
	rp30s15.shrexante.Q[[i]] <- unlist(randport.eval(rp30s15.Q[[i]], 
		additional.args=list(prices=us.pricemat[quarter.ends[i], ], 
		variance=varshr.Q[[i]]), keep="var.values"), use.names=FALSE)
}

# do graphs comparing ex ante factor model to ex ante shrinkage
P.facshr07Q4(NULL)
P.facshr08Q4(NULL)


# compute realized volatility for the quarters

rp30s15.realized.Q <- empty8
for(i in 1:8) {
	rp30s15.realized.Q[[i]] <- pp.realvol(valuation(rp30s15.Q[[i]], 
		us.pricemat[ quarter.ends[i]:quarter.ends[i+1], ], 
		collapse=TRUE))
}

# compute correlations by quarter of shrinkage versus realized
corshrreal <- numeric(8)
names(corshrreal) <- names(rp30s15.realized.Q)
corfacreal <- corshrreal

for(i in 1:8) {
	# ex ante is variance not volatility, 
	# but for Spearman only rank matters
	corshrreal[i] <- cor(rp30s15.shrexante.Q[[i]], 
		rp30s15.realized.Q[[i]], method='spearman')
	corfacreal[i] <- cor(rp30s15.facexante.Q[[i]], 
		rp30s15.realized.Q[[i]], method='spearman')
}

# do bootstrapping to find 95% confidence intervals for correlation

corboot.Q <- rep(list(numeric(5000)), 8)
names(corboot.Q) <- names(empty8)

for(i in 1:5000) {
	this.bootsamp <- sample(1000, 1000, replace=TRUE)
	for(j in 1:8) {
		corboot.Q[[j]][i] <- cor(rp30s15.shrexante.Q[[j]][
			this.bootsamp], rp30s15.realized.Q[[j]][
			this.bootsamp], method='spearman')
	}
}
corshrreal.95ci <- t(sapply(corboot.Q, quantile, probs=c(.025, .975)))

# do plot of correlations
P.corshrreal(NULL)


# dollar neutral portfolios

# create random portfolios
rpdn30s15.Q <- empty8
for(i in 1:8) {
	rpdn30s15.Q[[i]] <- random.portfolio(1000, 
		us.pricemat[quarter.ends[i], ], gross.value=1e6, net.value=0, 
		port.size=30, max.weight=.05, sum.weight=c("4"=.15))
}

# compute ex ante shrinkage volatility
rpdn30s15.shrexante.Q <- empty8
rpdn30s15.realized.Q <- rpdn30s15.shrexante.Q

for(i in 1:8) {
	rpdn30s15.shrexante.Q[[i]] <- unlist(randport.eval(rpdn30s15.Q[[i]], 
		additional.args=list(prices=us.pricemat[quarter.ends[i], ], 
		variance=varshr.Q[[i]]), keep="var.values"), use.names=FALSE)
}

# find realized 
for(i in 1:8) {
	# notice 'type' arguemnt to 'valuation'
	rpdn30s15.realized.Q[[i]] <- pp.realvol(valuation(rpdn30s15.Q[[i]], 
		us.pricemat[ quarter.ends[i]:quarter.ends[i+1], ], 
		collapse=TRUE, type="nav"))
}

# calculate quarterly correlations
corshrrealdn30 <- corshrreal
corshrrealdn30[] <- NA
for(i in 1:8) {
	corshrrealdn30[i] <- cor(rpdn30s15.shrexante.Q[[i]], 
		rpdn30s15.realized.Q[[i]], method='spearman')
}

# do plot of correlations
P.corshrrealdn(NULL)


# 120/20 portfolios

# create random portfolios
rp12020ls30s15.Q <- empty8
for(i in 1:8) {
	rp12020ls30s15.Q[[i]] <- random.portfolio(1000, 
		us.pricemat[quarter.ends[i], ], long.value=1.2e6, 
		short.value=2e5, port.size=30, max.weight=.05, 
		sum.weight=c("4"=.15))
}

# compute ex ante shrinkage volatility
rp12020ls30s15.shrexante.Q <- empty8
for(i in 1:8) {
	rp12020ls30s15.shrexante.Q[[i]] <- 
		unlist(randport.eval(rp12020ls30s15.Q[[i]], 
		additional.args=list(prices=us.pricemat[quarter.ends[i], ], 
		variance=varshr.Q[[i]]), keep="var.values"), use.names=FALSE)
}

# calculate realized
rp12020ls30s15.realized.Q <- empty8
for(i in 1:8) {
	# notice 'type' and 'cash' arguments to 'valuation'
	rp12020ls30s15.realized.Q[[i]] <- 
		pp.realvol(valuation(rp12020ls30s15.Q[[i]], 
		us.pricemat[ quarter.ends[i]:quarter.ends[i+1], ], 
		collapse=TRUE, type="nav", cash=0))
}

# calculate quarterly correlations
corshrreal12020ls <- corshrrealdn30
corshrreal12020ls[] <- NA
for(i in 1:8) {
	corshrreal12020ls[i] <- cor(rp12020ls30s15.shrexante.Q[[i]], 
		rp12020ls30s15.realized.Q[[i]], method='spearman')
}

# do plot of correlations
P.corshrreal12020ls(NULL)