Stock analysis and calculate VaR using R

This is a demo script for stock download and analysis using R

 (1) Download time series data from Yahoo Finance (example here is to download from 2011 up to 2013), or copy the url and download in browser and save as csv files

getstockdata.sh    Select all

curl "http://ichart.finance.yahoo.com/table.csv?s=MSFT&a=0&b=01&c=2011&d=11&e=31&f=2013&g=d" > msft.csv curl "http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=0&b=01&c=2011&d=11&e=31&f=2013&g=d" > aapl.csv curl "http://ichart.finance.yahoo.com/table.csv?s=^IXIC&a=0&b=01&c=2011&d=11&e=31&f=2013&g=d" > nasdaq.csv curl "http://ichart.finance.yahoo.com/table.csv?s=GOOG&a=0&b=01&c=2011&d=11&e=31&f=2013&g=d" > goog.csv ### or use wget "http://ichart.finance.yahoo.com/table.csv?s=MSFT&a=0&b=01&c=2011&d=11&e=31&f=2013&g=d" -O msft.csv

 (2) Run script in R, you can use source("stocks.r")

stocks.r    Select all

msft <- read.csv("msft.csv") aapl <- read.csv("aapl.csv") goog <- read.csv("goog.csv") nasdaq <- read.csv("nasdaq.csv") stocks <- nasdaq[,c("Date","Close")] colnames(stocks)[2] <- "NASDAQ" stocks["MSFT"] <- NA stocks$MSFT <- msft$Close stocks["AAPL"] <- NA stocks$AAPL <- aapl$Close stocks["GOOG"] <- NA stocks$GOOG <- goog$Close

 (3) Do some analysis in R

analysis.r    Select all

### plot Data plot(aapl$Date, aapl$Close) cor(goog[c("Open")],aapl[c("Open")]) summary(stocks) str(stocks) pairs(stocks[-1]) cor(stocks[-1]) cov(stocks[-1]) eigen(cor(stocks[-1])) ### Do linear regression mod1 <- lm(NASDAQ ~ MSFT + AAPL + GOOG, data=stocks) summary(mod1) ?summary.lm ### Confidence Limits on the Estimated Coefficients confint(mod1) summary.aov(mod1) ### Prediction of mean response for cases like this... predict(mod1, list(MSFT=30,AAPL=400,GOOG=800), interval="conf") ### Prediction for a single new case... predict(mod1, list(MSFT=30,AAPL=400,GOOG=800), interval="pred") mod2 = update(mod1,.~.-AAPL) summary(mod2) (prediction <- c(1,30,800) * coef(mod2)) sum(prediction) ### Compare two models anova(mod1, mod2) ### Regression Diagnostics par(mfrow=c(2,2)); plot(mod1); par(mfrow=c(1,1))

 (4) Stepwise Regression in R

stepwise.r    Select all

###retrieve more stocks data directly using R as below intc <- read.csv("http://ichart.finance.yahoo.com/table.csv?s=INTC&a=0&b=01&c=2011&d=11&e=31&f=2013&g=d") amzn <- read.csv("http://ichart.finance.yahoo.com/table.csv?s=AMZN&a=0&b=01&c=2011&d=11&e=31&f=2013&g=d") ###update model and do stepwise regression stocks$INTC <- intc$Close stocks$AMZN <- amzn$Close mod3 <- lm(NASDAQ ~ MSFT + AAPL + GOOG + INTC + AMZN, data=stocks) step(mod3, direction="both") summary(mod3) summary(step(mod3, direction="both"))

 (5) Delta Normal VaR and ES for a single asset using R
Reference (VaR using Excel) : http://www.youtube.com/watch?v=ZrKmVC-Ede8

DNVaR1.r    Select all

# get the time series of AAPL with Closed Price aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=2&b=23&c=2012&d=2&e=23&f=2013&g=d"); # calculate log returns, note the -log sign below aapl.logreturn = c(diff(-log(aapl$Close))); # calculate VaR for 95% confidence interval, VaR is 0.03513197 qnorm(0.95)*sd(aapl.logreturn)-mean(aapl.logreturn); #alternatively, qnorm(p=0.95, sd=sd(aapl.logreturn), mean=mean(aapl.logreturn)); # with a single asset of say $1M, the calculated Delta Normal VaR is $35,131.97 sprintf("%5.2f", 1e+06*(qnorm(0.95)*sd(aapl.logreturn)-mean(aapl.logreturn))); #Expected Shortfall(ES) for Normal Distribution = $62,227.24 v <- qnorm(p=0.95, sd=sd(aapl.logreturn), mean=mean(aapl.logreturn)) tailExp <- integrate(function(x) x * dnorm(x, mean=mean(aapl.logreturn), sd=sd(aapl.logreturn)), lower=-Inf, upper=v)$value / (1-0.95) -tailExp * 1e+06

 (6) Delta Normal VaR for a portfolio of assets using R
Reference (Portfolio VaRs using Excel) : https://www.youtube.com/watch?v=GqczSHRUaDk

DNVaR2.r    Select all

# get the time series of MSFT, GOOG, AAPL, ORCL and IBM msft <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=MSFT&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") goog <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=GOOG&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") orcl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=ORCL&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") ibm <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=IBM&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") # calculate log returns from Closed Prices msft.logreturn = c(diff(-log(msft$Close))); goog.logreturn = c(diff(-log(goog$Close))); aapl.logreturn = c(diff(-log(aapl$Close))); orcl.logreturn = c(diff(-log(orcl$Close))); ibm.logreturn = c(diff(-log(ibm$Close))); # calculate the individual volatility logreturns <- data.frame(msft.logreturn, goog.logreturn, aapl.logreturn, orcl.logreturn, ibm.logreturn); sapply(logreturns,sd); # calculate the individual VaR using Matrix and with position of $1M each using correlation matrix VaR1 = (1e+06*sapply(logreturns,sd)) %*% cor(logreturns); # calculate the Dollar Variance of the Portfolio VaR2 = VaR1 %*% (1e+06*sapply(logreturns,sd)); # Square Root Portfolio Variance to calculate the Delta Normal 1-day Portfolio VaR for 95% confidence interval. The result is $119,605.5 qnorm(0.95) * sqrt(VaR2); #Alternatively, the variance-cov method is demo here v1 = cov(logreturns) %*% rep(1e+06,5); v2 = rep(1e+06,5) %*% v1; sqrt(v2)*qnorm(0.95); # Portfolio VaR for 10 days sqrt(v2)*qnorm(0.95)* sqrt(10);

 (6.1) Delta Normal VaR and ES for a portfolio of assets using R
Using : simple return instead of log return

DNVaR2.r    Select all

# get the time series of MSFT, GOOG, AAPL, ORCL and IBM msft <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=MSFT&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") yhoo <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=YHOO&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") orcl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=ORCL&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") ibm <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=IBM&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") # calculate simple returns from Closed Prices #msft.simplereturn = -diff(msft$Close)/msft$Close[-length(msft$Close)]; #yhoo.simplereturn = -diff(yhoo$Close)/goog$Close[-length(yhoo$Close)]; #aapl.simplereturn = -diff(aapl$Close)/aapl$Close[-length(aapl$Close)]; #orcl.simplereturn = -diff(orcl$Close)/orcl$Close[-length(orcl$Close)]; #ibm.simplereturn = -diff(ibm$Close)/ibm$Close[-length(ibm$Close)]; msft.simplereturn = msft$Close[-nrow(msft)] / msft$Close[-1] -1; yhoo.simplereturn = yhoo$Close[-nrow(yhoo)] / yhoo$Close[-1] -1; aapl.simplereturn = aapl$Close[-nrow(aapl)] / aapl$Close[-1] -1; orcl.simplereturn = orcl$Close[-nrow(orcl)] / orcl$Close[-1] -1; ibm.simplereturn = ibm$Close[-nrow(ibm)] / ibm$Close[-1] -1; # plot the histogram, assume different stock position simplereturns <- data.frame(msft.simplereturn,yhoo.simplereturn,aapl.simplereturn,orcl.simplereturn,ibm.simplereturn); hist(as.matrix(simplereturns) %*% c(1e+06, 1e+05, 5e+03, 9e+04, 3e+05),col="green"); # calculate the individual volatility simplereturns <- data.frame(msft.simplereturn, yhoo.simplereturn, aapl.simplereturn, orcl.simplereturn, ibm.simplereturn); sapply(simplereturns,sd); #using the variance-covariance method with different stock position to calculate VaR $35,473.23 v1 = cov(simplereturns) %*% c(1e+06, 1e+05, 5e+03, 9e+04, 3e+05); v2 = c(1e+06, 1e+05, 5e+03, 9e+04, 3e+05) %*% v1; sqrt(v2)*qnorm(0.95); # Portfolio VaR for 10 days sqrt(v2)*qnorm(0.95)* sqrt(10); #calculate the ES $21,986.07 m <- sapply(simplereturns,mean) %*% c(1e+06, 1e+05, 5e+03, 9e+04, 3e+05); v <- qnorm(p=0.95, sd=sqrt(v2), mean=m); - integrate(function(x) x * dnorm(x, mean=m, sd=sqrt(v2)), lower=-Inf, upper=v)$value / (1-0.95) # use PerformanceAnalytics to calculate VaR and ES require(PerformanceAnalytics) stockportfolio.ts <- as.xts(data.frame(msft$Close,yhoo$Close,aapl$Close,orcl$Close,ibm$Close), order.by = as.Date(msft$Date,"%Y-%m-%d")) # calculate VaR $34,289.08 VaR(R = na.omit(Return.calculate(stockportfolio.ts, method="discrete")), p = 0.95, method = "gaussian", portfolio_method = "component", weights = c(1e+06, 1e+05, 5e+03, 9e+04, 3e+05)) # calculate ES $43,300.71 ETL(R = na.omit(Return.calculate(stockportfolio.ts, method="discrete")), p = 0.95, method = "gaussian", portfolio_method = "component", weights = c(1e+06, 1e+05, 5e+03, 9e+04, 3e+05))

 (7) Historical Simulation VaR for a single asset using R
Reference (Historical Simulation using Excel) : http://www.youtube.com/watch?v=6Nolb4-iRSI

HSVaR1.r    Select all

# get the time series of GOOG with Closed Price goog <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=GOOG&a=1&b=25&c=2008&d=6&e=18&f=2008&g=d") # calculate log returns, using natural log goog.logreturn = c(diff(-log(goog$Close))); # plot histogram hist(goog.logreturn, breaks=seq(-.11,.19,by=.01), col="green"); #show the first 10 items in ascending order head(sort(goog.logreturn),n=10); # find volatility for 5% percentile, that is -4.10867% quantile(goog.logreturn,0.05); #Historical Simulation VaR for $1,000 is 4.10867 -quantile(goog.logreturn,0.05)*1000; # find historical volatility for the 102 observations sd(goog.logreturn)*sqrt(length(goog[,1]));

 (7.1) Historical Simulation VaR for a single asset using R
Reference : Using the same data as per (5) above

HSVaR11.r    Select all

# get the time series of AAPL with Closed Price aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=2&b=23&c=2012&d=2&e=23&f=2013&g=d"); # calculate log returns, note the -log sign below aapl.logreturn = c(diff(-log(aapl$Close))); # calculate simple returns aapl.simplereturn = aapl$Close[-nrow(aapl)] / aapl$Close[-1] -1; # plot histogram and density curve hist(aapl.logreturn,breaks=seq(-.15,.1,by=.02), col="green",freq=FALSE,main="Asset Returns", xlab="Return%",xlim=range(-0.15,0.1),ylim=range(0,28));lines(density(aapl.logreturn),col="red", lty="dotted",panel.last=abline(v=quantile(aapl.logreturn, 0.05), col="darkred",lwd=3)); mtext("95% C.I.",side=3,line=-3,adj=0.40,col="darkred"); #show the first 15 items in ascending order head(sort(aapl.logreturn),n=15); # find VaR for 5% percentile, volatility is 0.03205498 = $32,054.98 -quantile(aapl.logreturn,0.05)*1e+06; # find VaR for 5% percentile, volatility is 0.03205498 = $31,546.66 -quantile(aapl.simplereturn,0.05)*1e+06; # with a single asset of say $1M, the calculated Historical Simulation VaR is $32,054.98 comparing with the result of Delta Normal VaR of $35,131.97 as per (5) above sprintf("%5.2f", 1e+06*-quantile(aapl.logreturn,0.05)); # calculate ES $46,965.9 for log returns -mean(aapl.logreturn[aapl.logreturn <= quantile(aapl.logreturn, 0.05)])*1e+6; # calculate ES $45,565.73 for simple returns -mean(aapl.simplereturn[aapl.simplereturn <= quantile(aapl.simplereturn, 0.05)])*1e+6;

 (7.2) Historical Simulation VaR for a single asset using R and Package PerformanceAnalytics
Reference : http://cran.at.r-project.org/web/packages/PerformanceAnalytics/
http://cran.r-project.org/web/packages/quantmod/

HSVaR12.r    Select all

#install.packages("PerformanceAnalytics"); require(PerformanceAnalytics) # get the time series of AAPL with Closed Price aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=2&b=23&c=2012&d=2&e=23&f=2013&g=d"); # convert to time series aapl.ts <- as.xts(aapl$Close, order.by = as.Date(aapl$Date,"%Y-%m-%d")); # calculate VaR $32,054.98 for log returns -VaR(Return.calculate(aapl.ts, method="log"),p=.95,method="historical") * 1e+06; # calculate VaR $31,546.66 for simple returns -VaR(Return.calculate(aapl.ts, method="discrete"),p=.95,method="historical") * 1e+06; # calculate ES $46,965.9 for log returns -ETL(Return.calculate(aapl.ts, method="log"),p=.95,method="historical") * 1e+06; # calculate ES $45,565.73 for simple returns -ETL(Return.calculate(aapl.ts, method="discrete"),p=.95,method="historical") * 1e+06; #plot PerformanceSummary using Package PerformanceAnalytics charts.PerformanceSummary(Return.calculate(aapl.ts, method="log")); #plot the time series chart using Package quantmod #install.packages("quantmod"); require("quant mod"); chartSeries(aapl.ts, theme = chartTheme("white"), TA = c(addBBands(),addTA(RSI(aapl.ts))));

 (7.3) Calculate VaR for a single asset using Monte carlo simulation: Brownian motion
Require : http://cran.at.r-project.org/web/packages/PerformanceAnalytics/

MCVaR.r    Select all

#install.packages("PerformanceAnalytics"); require(PerformanceAnalytics) # get the time series of AAPL with Closed Price aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=2&b=23&c=2012&d=2&e=23&f=2013&g=d"); # convert to time series aapl.ts <- as.xts(aapl$Close, order.by = as.Date(aapl$Date,"%Y-%m-%d")); # calculate VaR $32,054.98 for log returns -VaR(Return.calculate(aapl.ts, method="log"),p=.95,method="historical") * 1e+06; # calculate VaR $3x,xxx (will print 3 iterations) using random 250 samples log returns generated based on sd and mean of AAPL aapl.logreturn = c(diff(-log(aapl$Close))); for (i in 1:3) {print (-VaR(as.xts(rnorm(250, sd=sd(aapl.logreturn), mean=mean(aapl.logreturn)),order.by=as.Date(aapl$Date,"%Y-%m-%d")),p=.95,method="historical") * 1e+06)}; # calculate VaR using random 1000 samples log returns generated based on sd and mean of AAPL for (i in 1:3) {print (-VaR(as.xts(rnorm(1000, sd=sd(aapl.logreturn), mean=mean(aapl.logreturn)),order.by=as.Date(Sys.Date():(Sys.Date()+1000-1))),p=.95,method="historical") * 1e+06)}; # calculate VaR (will print 3 iterations) using random 250 samples log returns generated # based on the assumptions of annual return drift=0.10 and annual return volatility=0.40 # daily volatility = annual volatility / sort(T) # daily mean drift = annual drift/T - 0.5 * sd^2 # model reference from http://www.youtube.com/watch?v=e79OtCamxD0 for (i in 1:3) {print (-VaR(as.xts(rnorm(250, sd=0.40/sqrt(250), mean=(0.10/250)-0.5*((0.40/sqrt(250))^2)),order.by=as.Date(aapl$Date,"%Y-%m-%d")),p=.95,method="historical") * 1e+06)}; # calculate VaR (will print 3 iterations) using random 1000 samples log returns generated # based on the assumptions of annual return drift=0.10 and a higher annual return volatility=0.60 for (i in 1:3) {print (-VaR(as.xts(rnorm(1000, sd=0.60/sqrt(250), mean=(0.10/250)-0.5*((0.60/sqrt(250))^2)),order.by=as.Date(Sys.Date():(Sys.Date()+1000-1))),p=.95,method="historical") * 1e+06)}; # plot the simulated prices from one instance of simulation aapl.mc.logreturns = rnorm(1000, sd=0.60/sqrt(250), mean=(0.10/250)-0.5*((0.60/sqrt(250))^2)); plot(as.xts(c(1,exp(cumsum(aapl.mc.logreturns))) * tail(aapl$Close,n=1),order.by=as.Date(as.Date(c("2010-01-01")):(as.Date(c("2010-01-01"))+1000))),main="Simulated Price");

 (8) Historical Simulation VaR for a portfolio of assets using R
Reference : Using the same data as per (6) Delta Normal VaR for a portfolio of assets

HSVaR2.r    Select all

# get the time series of MSFT, GOOG, AAPL, ORCL and IBM msft <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=MSFT&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") goog <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=GOOG&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") orcl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=ORCL&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") ibm <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=IBM&a=5&b=14&c=2011&d=5&e=8&f=2012&g=d") # calculate log returns from Closed Prices msft.logreturn = c(diff(-log(msft$Close))); goog.logreturn = c(diff(-log(goog$Close))); aapl.logreturn = c(diff(-log(aapl$Close))); orcl.logreturn = c(diff(-log(orcl$Close))); ibm.logreturn = c(diff(-log(ibm$Close))); # plot the histogram and density logreturns <- data.frame(msft.logreturn, goog.logreturn, aapl.logreturn, orcl.logreturn, ibm.logreturn); hist(as.matrix(logreturns) %*% rep(1e+06,5),col="green",freq=FALSE); lines(density(as.matrix(logreturns) %*% rep(1e+06,5)),col="red") # find VaR for 5% percentile, that is $109,409.9 and comparing with the result of Delta Normal VaR of $119,605.5 as per (6) above -quantile(as.matrix(logreturns) %*% rep(1e+6,5), 0.05);

 (8.1) Historical Simulation VaR and ES for a portfolio of assets using R
Using : simple return instead of log return
Calculate : Incremental VaR, Marginal VaR and Component VaR

HSVaR250.r    Select all

# get the time series of MSFT, YHOO, AAPL, ORCL and IBM for 251 trading days msft <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=MSFT&a=3&b=3&c=2013&d=2&e=31&f=2014&g=d") yhoo <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=YHOO&a=3&b=3&c=2013&d=2&e=31&f=2014&g=d") aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=3&b=3&c=2013&d=2&e=31&f=2014&g=d") orcl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=ORCL&a=3&b=3&c=2013&d=2&e=31&f=2014&g=d") ibm <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=IBM&a=3&b=3&c=2013&d=2&e=31&f=2014&g=d") # calculate simple returns from Closed Prices msft.simplereturn = msft$Close[-nrow(msft)] / msft$Close[-1] -1; yhoo.simplereturn = yhoo$Close[-nrow(yhoo)] / yhoo$Close[-1] -1; aapl.simplereturn = aapl$Close[-nrow(aapl)] / aapl$Close[-1] -1; orcl.simplereturn = orcl$Close[-nrow(orcl)] / orcl$Close[-1] -1; ibm.simplereturn = ibm$Close[-nrow(ibm)] / ibm$Close[-1] -1; # plot the histogram simplereturns <- data.frame(msft.simplereturn,yhoo.simplereturn,aapl.simplereturn,orcl.simplereturn,ibm.simplereturn); hist(as.matrix(simplereturns) %*% rep(1,5),col="green",freq=FALSE); lines(density(as.matrix(simplereturns) %*% rep(1,5)),col="red"); # find VaR for 5% percentile, that is $71,726.65 -quantile(as.matrix(simplereturns) %*% rep(1e+6,5), 0.05); # find Incremental VaR for each asset by removing each asset from the portfolio # Incremental VaR : MSFT $3396.407; YHOO $19522.17; AAPL $10691.23; ORCL $12628.08; IBM $10749.8 Iport = matrix(1e+06, nrow=5, ncol=5); diag(Iport)=0; for (i in 1:5) {print(quantile(as.matrix(simplereturns) %*% (Iport)[i,], 0.05) -quantile(as.matrix(simplereturns) %*% rep(1e+6,5), 0.05))}; # Marginal VaR by list out the 13th element sort by portfolio return column # msft.simplereturn yhoo.simplereturn aapl.simplereturn orcl.simplereturn ibm.simplereturn portfolio #45 -0.02118989 -0.03323661 0.008112513 -0.01670709 -0.009686039 -0.07270712 returns <- data.frame(simplereturns, portfolio=as.matrix(simplereturns) %*% rep(1,5) ) tail(head(returns[order(returns$portfolio),],n=length(returns[,1])*0.05+1),n=1) # Component VaR by multiplying Marginal VaR by the weight of each asset # msft.simplereturn yhoo.simplereturn aapl.simplereturn orcl.simplereturn ibm.simplereturn portfolio #45 -21189.89 -33236.61 8112.513 -16707.09 -9686.039 -72707.12 tail(head(returns[order(returns$portfolio),],n=250*0.05+1),n=1) * 1e+6 #find Expected Shortfall (ES) that is $103,035.5 m <- as.matrix(simplereturns) %*% rep(1e+6,5); v <- quantile(m, 0.05); -mean(m[m <= v]); # plot all hist(as.matrix(simplereturns) %*% rep(1,5),col="green",freq=FALSE,main="Portfolio Returns", xlab="Return%",xlim=range(-0.2,0.15),ylim=range(0,10));lines(density(as.matrix(simplereturns) %*% rep(1,5)),col="red", lty="dotted",panel.last=abline(v=quantile(as.matrix(simplereturns) %*% rep(1,5), 0.05), col="darkred",lwd=3)); mtext("95% C.I.",side=3,line=-3,adj=0.28,col="darkred"); #plot all plus superimposed normal distribution curve xvals <- seq(from=-0.15, to=0.15,length=100) hist(as.matrix(simplereturns) %*% rep(1,5),col="green",freq=FALSE,main="Portfolio Returns", xlab="Return%",xlim=range(-0.2,0.15),ylim=range(0,10));lines(density(as.matrix(simplereturns) %*% rep(1,5)),col="red", lty="dotted",panel.last=abline(v=quantile(as.matrix(simplereturns) %*% rep(1,5), 0.05), col="darkred",lwd=3)); mtext("95% C.I.",side=3,line=-3,adj=0.28,col="darkred");lines(xvals,dnorm(xvals,mean(as.matrix(simplereturns) %*% rep(1,5)),sd(as.vector(as.matrix(simplereturns) %*% rep(1,5)))),lwd=1,col="blue");

 (8.2) Historical Simulation 10-Days VaR and ES for a portfolio of assets using R
Using : log return for 10 days and simple return for portfolio with 4 years historical data

HSVaR10.r    Select all

# get the time series of MSFT, YHOO, AAPL, ORCL and IBM for 260 trading days msft <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=MSFT&a=2&b=20&c=2010&d=2&e=31&f=2014&g=d") yhoo <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=YHOO&a=2&b=20&c=2010&d=2&e=31&f=2014&g=d") aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=2&b=20&c=2010&d=2&e=31&f=2014&g=d") orcl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=ORCL&a=2&b=20&c=2010&d=2&e=31&f=2014&g=d") ibm <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=IBM&a=2&b=20&c=2010&d=2&e=31&f=2014&g=d") # calculate log returns from Closed Prices and takes first 1000 samples only msft.logreturn10 = c(diff(-log(msft$Close),10)); yhoo.logreturn10 = c(diff(-log(yhoo$Close),10)); aapl.logreturn10 = c(diff(-log(aapl$Close),10)); orcl.logreturn10 = c(diff(-log(orcl$Close),10)); ibm.logreturn10 = c(diff(-log(ibm$Close),10)); logreturns10 <- data.frame(msft.logreturn10, yhoo.logreturn10, aapl.logreturn10, orcl.logreturn10, ibm.logreturn10)[1:1000,]; # convert to simple returns simplereturns10 <- (exp(logreturns10)-1) # find VaR for 5% percentile for 10 days holding period, that is $259,314.5 -quantile(as.matrix(simplereturns10) %*% rep(1e+6,5), 0.05); #find 10-Days expected Shortfall (ES) that is $343,451.5 m <- as.matrix(simplereturns10) %*% rep(1e+6,5); v <- quantile(m, 0.05); -mean(m[m <= v]); # plot all hist(as.matrix(simplereturns10) %*% rep(1,5),col="green",freq=FALSE,main="Portfolio Returns for Holding Period 10 days", xlab="Return%",xlim=range(-0.6,0.6),ylim=range(0,4));lines(density(as.matrix(simplereturns10) %*% rep(1,5)),col="red", lty="dotted",panel.last=abline(v=quantile(as.matrix(simplereturns10) %*% rep(1,5), 0.05), col="darkred",lwd=3)); mtext("95% C.I.",side=3,line=-3,adj=0.20,col="darkred"); #List the 51st item sort by portfolio return #587 0.004545455 0.00523416 -0.002571286 -0.003148254 -0.01536958 -0.2592768 returns10 <- data.frame(simplereturns, portfolio=as.matrix(simplereturns10) %*% rep(1,5)); tail(head(returns10[order(returns10$portfolio),],n=length(returns10[,1])*0.05+1),n=1);

 (8.3) Historical Simulation VaR and ES for a portfolio of assets using R and Package PerformanceAnalytics
Backtesting : if Actual Return < Theoretical VaR
constrOptim Ref : http://www.youtube.com/watch?v=MCvz-c6UUkw

HSVaR21.r    Select all

require(PerformanceAnalytics) # get the time series of MSFT, YHOO, AAPL, ORCL and IBM msft <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=MSFT&a=0&b=1&c=2012&d=2&e=31&f=2014&g=d") yhoo <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=YHOO&a=0&b=1&c=2012&d=2&e=31&f=2014&g=d") aapl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=AAPL&a=0&b=1&c=2012&d=2&e=31&f=2014&g=d") orcl <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=ORCL&a=0&b=1&c=2012&d=2&e=31&f=2014&g=d") ibm <-read.csv("http://ichart.finance.yahoo.com/table.csv?s=IBM&a=0&b=1&c=2012&d=2&e=31&f=2014&g=d") stockportfolio.ts <- as.xts(data.frame(msft$Close,yhoo$Close,aapl$Close,orcl$Close,ibm$Close), order.by = as.Date(msft$Date,"%Y-%m-%d")) # calculate VaR 0.0215341 0.01744393 0.02875233 0.02034777 0.01493424 -VaR(Return.calculate(stockportfolio.ts[1:251], method="discrete"),p=.95,method="historical") # calculate ES 0.02622294 0.02608046 0.03779529 0.02482448 0.02314859 -ETL(Return.calculate(stockportfolio.ts[1:251], method="discrete"),p=.95,method="historical") # Backtesting p=0.95 assuming equal weights # i Actual Theoretical # 2013-06-20 368 -0.128825 0.1031542 # 2013-06-21 369 -0.1229124 0.1039837 # 2013-07-19 388 -0.1754735 0.103187 # 2013-08-27 415 -0.1080645 0.1028595 # 2014-02-03 524 -0.1150287 0.1047768 for (i in 1:(nrow(stockportfolio.ts)-251)) { actual <- sum(na.omit(Return.calculate(stockportfolio.ts[(251+i-1):(251+i)]))); theoretical <- sum(-VaR(Return.calculate( stockportfolio.ts[(i):(251+i-1)],method="discrete"),p=.95,method="historical")); if (any(actual < -theoretical)) { print(cbind(stockportfolio.ts[(251+i),0],i=251+i,Actual=actual,Theoretical=theoretical)) } } # Backtesting p=0.99 assuming equal weights # i Actual Theoretical # 2013-07-19 388 -0.1754735 0.1667104 for (i in 1:(nrow(stockportfolio.ts)-251)) { actual <- sum(na.omit(Return.calculate(stockportfolio.ts[(251+i-1):(251+i)]))); theoretical <- sum(-VaR(Return.calculate( stockportfolio.ts[(i):(251+i-1)],method="discrete"),p=.99,method="historical")); if (any(actual < -theoretical)) { print(cbind(stockportfolio.ts[(251+i),0],i=251+i,Actual=actual,Theoretical=theoretical)) } } # Backtesting p=0.99 assuming weights of c(30,50,200,50,50) # i Actual Theoretical # 2013-01-24 266 -23.5823 13.22212 # 2014-01-28 520 -13.10302 12.38539 for (i in 1:(nrow(stockportfolio.ts)-251)) { actual <- sum(na.omit(Return.calculate(stockportfolio.ts[(251+i-1):(251+i)])) * c(30,50,200,50,50)); theoretical <- sum(-VaR(Return.calculate( stockportfolio.ts[(i):(251+i-1)],method="discrete"),p=.99,method="historical") * c(30,50,200,50,50)); if (any(actual < -theoretical)) { print(cbind(stockportfolio.ts[(251+i),0],i=251+i,Actual=actual,Theoretical=theoretical)) } } # Use constraint optimisation to find the weight of each share of the portfolio # optimise function g and to minimise VaR portfolioVaR <- -VaR(Return.calculate(stockportfolio.ts[1:251], method="discrete"),p=.95,method="historical"); g <- function(x) {sum(portfolioVaR* x)}; # constraint matrix A <- matrix(c(rep(1,5),rep(-1,5),diag(1,5),diag(-1,5)),nrow=2+5+5, ncol=5, byrow=TRUE); # constraint values : value of portfolio = 500 and shareholding for each stock > 0 that is no short selling b <- 500 * c(0.999999,-1.000-0.000001,rep(0.000001,5),rep(-1.000-0.000001,5)); # run constrOptim with initial value of 100 each # $par = 0.1113882 58.9182990 0.4609939 78.1495814 362.3602354 # $value = 8.045164 constrOptim(theta=rep(100,5),grad=NULL,f=g,ui=A,ci=b) # Backtesing again and plot the graph for Actual vs. Theoretical results based on the above optimised portfolio # i Actual Theoretical # 2013-04-19 325 -28.85769 12.51676 # 2013-10-17 451 -24.06681 13.6749 result.frame = data.frame() for (i in 1:(nrow(stockportfolio.ts)-251)) { actual <- sum(na.omit(Return.calculate(stockportfolio.ts[(251+i-1):(251+i)])) * c(0.1113882,58.9182990,0.4609939,78.1495814,362.3602354)); theoretical <- sum(-VaR(Return.calculate( stockportfolio.ts[(i):(251+i-1)],method="discrete"),p=.99,method="historical") * c(0.1113882,58.9182990,0.4609939,78.1495814,362.3602354)); if (any(actual < -theoretical )) { print(cbind(stockportfolio.ts[(251+i),0],i=251+i,Actual=actual,Theoretical=theoretical)); }; result.frame <- rbind(result.frame, data.frame(stockportfolio.ts[(251+i),],i=(251+i), Actual=actual,Theoretical=-theoretical)); } result.zoo <- as.zoo(result.frame) plot(x = result.zoo[,7:8], screens = 1, main="Actual vs. Theoretical VaR", xlab = "Time", ylab = "Return", col = c("darkgreen", "red")); legend(x = "bottomright", legend = c("Actual","Theoretical"),lty=1,col=c("darkgreen", "red")); # revise the optimisation function to include portfolioLogReturn and to minimise risk / return ratio #$par 1.177901 111.451867 142.695042 143.658599 101.017062 portfolioLogReturn <- colSums(Return.calculate(stockportfolio.ts[1:251], method="log"),na.rm=TRUE) g1 <- function(x) {sum(portfolioVaR * x)/sum(portfolioLogReturn * x)}; A <- matrix(c(rep(1,5),rep(-1,5),diag(1,5),diag(-1,5)),nrow=2+5+5, ncol=5, byrow=TRUE); b <- 500 * c(0.999999,-1.000-0.000001,rep(0.000001,5),rep(-1.000-0.000001,5)); constrOptim(theta=rep(100,5),grad=NULL,f=g1,ui=A,ci=b) # rerun backtesting with portfolio weights = c(1.177901,111.451867,142.695042,143.658599,101.017062) result.frame = data.frame() for (i in 1:(nrow(stockportfolio.ts)-251)) { actual <- sum(na.omit(Return.calculate(stockportfolio.ts[(251+i-1):(251+i)])) * c(1.177901,111.451867,142.695042,143.658599,101.017062)); theoretical <- sum(-VaR(Return.calculate( stockportfolio.ts[(i):(251+i-1)],method="discrete"),p=.99,method="historical") * c(1.177901,111.451867,142.695042,143.658599,101.017062)); if (any(actual < -theoretical )) { print(cbind(stockportfolio.ts[(251+i),0],i=251+i,Actual=actual,Theoretical=theoretical)); }; result.frame <- rbind(result.frame, data.frame(stockportfolio.ts[(251+i),],i=(251+i), Actual=actual,Theoretical=-theoretical)); } result2.zoo <- as.zoo(result.frame) plot(x = result2.zoo[,7:8], screens = 1, main="Actual vs. Theoretical VaR with optimised portfolio", xlab = "Time", ylab = "Return", col = c("darkgreen", "red")); legend(x = "topleft", legend = c("Actual","Theoretical"),lty=1,col=c("darkgreen", "red"));

There is one important difference between log return and simple return is that log return is time additive and simple return is portfolio additive. (Reference : http://www.youtube.com/watch?v=PtoUlt3V0CI)

How to upgrade svn to 1.7 for Mac OS X Mountain Lion

upgradesvn.sh    Select all

svn --version
cd ~/Downloads/
curl -o subversion-latest.tar.gz http://apache.mirrors.tds.net/subversion/subversion-1.7.11.tar.gz
tar -xvf subversion-latest.tar.gz
cd subversion-1.7.11/
sh get-deps.sh neon
cd neon/
./configure --with-ssl
make
sudo make install
cd ..
./configure --prefix=/usr --with-neon
make
sudo make install
svn --version


One of the major changes in this release 1.7 is a move to a new metadata format that does not require multiple .svn directories in working copies. Instead, working copies will have a single metadata directory.

How to install watir on Windows XP / 7 (32 bit) for IE Webapp testing

(1) Goto http://rubyinstaller.org/downloads and download
Ruby 1.9.3-p429
DevKit-tdm-32-4.5.2-20111229-1559-sfx.exe


(2) Download IEDriverServer_Win32_2.33.0.zip from
http://code.google.com/p/selenium/downloads/list
and extract to C:\Windows\


(3) install Ruby 1.9.3-p429 using Administrator
Choose Add Ruby executables to PATH and Associate .rb and .rbw
Install to say C:\Ruby193

(4) Double click and extract DevKit-tdm to say C:\Ruby193\DevKit


(5) Use Command Prompt (cmd.exe) and enter into DOS shell and run
cd C:\Ruby193\DevKit
ruby dk.rb init
ruby dk.rb install
gem install watir --no-ri --no-rdoc
# might have some errors after the last step
gem install nokogiri
gem install mini_magick -v 3.5
gem install watir --no-ri --no-rdoc
gem install watir-webdriver --no-ri --no-rdoc


(6) Use irb to test interative session on IE
require "watir"
browser = Watir::Browser.new
browser.goto("www.google.com")
browser.text_field(:name => 'q').set 'Watir Example'
browser.button(:name => 'btnK').click
browser.close

require "watir-webdriver"
browser = Watir::Browser.new :ie
browser.goto 'http://bit.ly/watir-example'
browser.text_field(:name => 'entry.0.single').set 'Watir'
browser.radio(:value => 'Watir').set

Personal Migration Guide for Google App Engine from python 2.5 to 2.7

app.yaml    Select all

#change application name to new id that enable HRD #application: myapp application: myapp-hrd #change runtime #runtime: python runtime: python27 #add threadsafe threadsafe: true #change remote_api handlers: #- url: /remote_api # script: $PYTHON_LIB/google/appengine/ext/remote_api/handler.py # login: admin - url: /remoteapi.* script: google.appengine.ext.remote_api.handler.app login: admin #change url for app - url: .* script: myapp.app # script: myapp.py #add libraries libraries: - name: webapp2 version: "2.5.2" - name: jinja2 version: "2.6" #add services inbound_services: - warmup

myapp.py    Select all

#import webapp2 and jinja2 import webapp2 import jinja2 #import wsgiref.handlers #from google.appengine.ext import webapp #from google.appengine.ext.webapp import RequestHandler #from google.appengine.ext.webapp import template #from google.appengine.ext.webapp.util import run_wsgi_app #use ndb instead of db #from google.appengine.ext import db #class Greeting(db.Model): # author = db.UserProperty() # content = db.StringProperty(indexed=False) # role = db.StringProperty(required=True, choices=set(["executive", "manager", "helper"])) # date = db.DateTimeProperty(auto_now_add=True) from google.appengine.ext import ndb class Greeting(ndb.Model): """Models an individual Guestbook entry with author, content, and date.""" author = ndb.UserProperty() content = ndb.StringProperty(indexed=False) role = ndb.StringProperty(required=True, choices=set(["executive", "manager", "helper"])) date = ndb.DateTimeProperty(auto_now_add=True) #use TEMPLATE_DIR and JINJA_ENVIRONMENT TEMPLATE_DIR = 'templates' JINJA_ENVIRONMENT = jinja2.Environment( loader=jinja2.FileSystemLoader(os.path.join(os.path.dirname(__file__), TEMPLATE_DIR)), # loader=jinja2.FileSystemLoader(os.path.dirname(__file__)), extensions=['jinja2.ext.autoescape']) #use jinja2 template # path = os.path.abspath(os.path.join(os.path.dirname(__file__), TEMPLATE_DIR, fileName)) # self.response.out.write(template.render(path, values)) template = JINJA_ENVIRONMENT.get_template(fileName) self.response.write(template.render(values)) #use webapp2 instead of webapp #def main(): #application = webapp.WSGIApplication( app = webapp2.WSGIApplication( [('/', MainHandler) # %3D does not work now # ,(r'/id%3D([A-Za-z0-9-.]*)', IDHandler) ,(r'/id=([A-Za-z0-9-.]*)', IDHandler) ], debug=True) # run_wsgi_app(application) #if __name__ == '__main__': # main()

Then, Duplicate Application Settings to a new application id say myapp-hrd

Finally, Migrate the app to the High Replication Datastore (HRD) in Google App Engine admin before uploading the python 2.7 code to the newly converted application id myapp-hrd

At the end of the migration process, the new HRD app (myapp-hrd) will be aliased to the old appid.

How to enable Java 6 Web plug-in for Max OS X Lion / Mountain Lion / Mavericks

http://support.apple.com/kb/ht5559

JAVA SE6 (OS X v10.7/10.8/10.9) -> http://support.apple.com/kb/DL1572
JAVA SE6 (OS X v10.6) -> http://support.apple.com/kb/DL1573
JAVA SE5 & SE6 (OS X v10.5) -> http://support.apple.com/kb/DL1359