## Thursday, 31 January 2019

### Polynomials-R

Polynomials-R
> install.packages("ar.matrix")
require("ggplot2")
> obs <- r.AR1(100, M=30, sigma=1, rho=.98)
> dim(obs)
 100  30
> obs_adj <- obs - obs[,1]
> ar1_df <- data.frame(obs=c(t(obs_adj)), realization=rep(1:100, each=30),
+                      time=rep(1:30, 100))
> ggplot(data=ar1_df, aes(time, obs, group=realization, color=realization)) +
+     geom_line()

> install.packages("mvp")
coeffs(a) <- 1
> coeffs(a) <- 0
> p <- rmvp(10,9,9,letters[1:4])
> deriv(p,letters[1:3])
mvp object algebraically equal to
2448 a^3 b^3 c^16 d^17  +  2145 a^4 b^12 c^10 d^7  +  4400 a^10 b^19 c d^18  +  14280 a^16 b^14 c^7 d^15  +  11628 a^18 b^5 c^16 d^4  +  2025 a^26 b^4 c^2 d
> deriv(p,rev(letters[1:3]))
> deriv(p,rev(letters[1:3]))
> x <- rmvp(7,symbols=6)
> v <- allvars(x)
> s <- as.mvp("1  +  y  -  y^2 zz  +  y^3 z^2")
> LHS <- subsmvp(deriv(x,v)*deriv(s,"y"),v,s)
> RHS <- deriv(subsmvp(x,v,s),"y")
> LHS - RHS
mvp object algebraically equal to
0
> p <- as.mvp("1+a^2 + a*b^2 + c")
> p
mvp object algebraically equal to
1  +  a b^2  +  a^2  +  c
> f <- as.function(p)
> f(a=1)
mvp object algebraically equal to
2  +  b^2  +  c
> f(a=1,b=2)
mvp object algebraically equal to
6  +  c
> f(a=1,b=2,c=3)
 9
> f(a=1,b=2,c=3,drop=FALSE)
mvp object algebraically equal to
9
> x == mpoly_to_mvp(mpoly::as.mpoly(x))
 TRUE
> kahle  <- mvp(
+     vars   = split(cbind(letters,letters[c(26,1:25)]),rep(seq_len(26),each=2)),
+     powers = rep(list(1:2),26),
+     coeffs = 1:26
+ )
library(mpoly)
> f<- mp("1-2x + x^2 + 100 x^2 y +100y^2")
> f <- as.function(f)
f(.) with . = (x, y)
> df <- expand.grid(x=seq(-2, 2,.01),y=seq(-1, 3, .01))
library(scales)
qplot(x , y, data = df,   geom = c("raster","contour"),fill = f+.001, z=f, colour= I("white"),bins=6) + scale_fill_gradient2(low = "blue",mid = "red",high = "green", trans="log10", guide = "none")

## Wednesday, 30 January 2019

### Equity Return Calculation

Equity Return Calculation
SieDates <- as.character(format(as.POSIXct(attr(siemens,"times")),"%Y-%m-%d"))
SieRet <- timeSeries(siemens * 100,charvec = SieDates)
colnames(SieRet) <- "SieRet"
## Stylised Facts 1
par(mfrow = c(2,2))
seriesPlot(SieRet,title = FALSE,main = "Daily Returns of Siemens",col = "blue")
boxPlot(SieRet, title = FALSE, main = "Box plot of Returns",col = "red",cex =0.5,pch = 19)
acf(SieRet,main ="ACF of Returns",lag.max = 20, ylab = "", xlab ="",col = "blue",ci.col = "red")
pacf(SieRet,main ="PACF of Returns",lag.max = 20,ylab ="",xlab= "",col="blue",ci.col="red")

## Stylised Facts II
SieRetAbs <- abs(SieRet)
SieRet100 <- tail(sort(abs(series(SieRet))), 100)
idx <- which(series(SieRetAbs) > SieRet100,arr.ind =TRUE)
SieRetAbs100 <- timeSeries(rep(0,length(SieRet)),charvec = time(SieRet))
SieRetAbs100[idx , 1] <- SieRetAbs[idx]
acf(SieRetAbs,main="ACF of Absolute Returns",lag.max = 20,ylab="",xlab="",col ="blue",ci.col="red")
pacf(SieRetAbs,main="PACF of Absolute Returns",lag.max = 20, ylab="", xlab="",col="blue",ci.col="red")
qqnormPlot(SieRet,main="QQ-Plot of Returns",title = FALSE,col = "blue",cex = 0.5,pch= 19)
plot(SieRetAbs100,type="h",main="Volatility Clustering",ylab="",xlab="",col="blue")

## Monday, 28 January 2019

### Portfolio

Portfolio
library("fPortfolio", lib.loc="~/R/win-library/3.5")
Warning messages:
1: package ‘fPortfolio’ was built under R version 3.5.2
2: package ‘timeDate’ was built under R version 3.5.1
3: package ‘timeSeries’ was built under R version 3.5.1
4: package ‘fBasics’ was built under R version 3.5.1
5: package ‘fAssets’ was built under R version 3.5.2
> library(fPortfolio)
> donlp2NLPControl(
+     iterma=4000, nstep=20, fnscale=1, report=FALSE, rep.freq=1,
+     tau0=1, tau=0.1, del0=1, epsx=1e-05, delmin=0.1 * del0,
+     epsdif=1e-08, nreset.multiplier=1, difftype=3, epsfcn=1e-16,
+     taubnd=1, hessian=FALSE, te0=TRUE, te1=FALSE, te2=FALSE,
+     te3=FALSE, silent=TRUE, intakt=TRUE)
Error in donlp2NLPControl(iterma = 4000, nstep = 20, fnscale = 1, report = FALSE,  :
> rdonlp2NLP(start, objective,
+            lower=0, upper=1, linCons, funCons, control=list())
Error in rdonlp2NLP(start, objective, lower = 0, upper = 1, linCons, funCons,  :
> donlp2NLP(start, objective,
+           par.lower=NULL, par.upper=NULL,
+           eqA=NULL, eqA.bound=NULL,
+           ineqA=NULL, ineqA.lower=NULL, ineqA.upper=NULL,
+           eqFun=list(), eqFun.bound=NULL,
+           ineqFun=list(), ineqFun.lower=NULL, ineqFun.upper=NULL,
+           control=list())
Error in loadNamespace(name) : there is no package called ‘Rdonlp2’
> nlminb2NLPControl(
+     eval.max=500, iter.max=400, trace=0, abs.tol=1e-20, rel.tol=1e-10,
+     x.tol=1.5e-08, step.min=2.2e-14, scale=1, R=1, beta.tol=1e-20)
\$`eval.max`
 500

\$iter.max
 400

\$trace
 0

\$abs.tol
 1e-20

\$rel.tol
 1e-10

\$x.tol
 1.5e-08

\$step.min
 2.2e-14

\$scale
 1

\$R
 1

\$beta.tol
 1e-20

> rnlminb2
function (...)
{
Rnlminb2::nlminb2(...)
}
<bytecode: 0x000000001d797450>
<environment: namespace:fPortfolio>
> ramplNLP(start, objective,
+          lower=0, upper=1, amplCons, control=list(), ...)
Error: '...' used in an incorrect context
> amplNLP()
 NA
> amplNLPControl(
+     solver="minos", project="ampl", trace=FALSE)
\$`solver`
 "minos"

\$project
 "ampl"

\$trace
 FALSE

> rdonlp2NLP(start, objective,
+            lower=0, upper=1, linCons, funCons, control=list())
Error in rdonlp2NLP(start, objective, lower = 0, upper = 1, linCons, funCons,  :
> donlp2NLP(start, objective,
+           par.lower=NULL, par.upper=NULL,
+           eqA=NULL, eqA.bound=NULL,
+           ineqA=NULL, ineqA.lower=NULL, ineqA.upper=NULL,
+           eqFun=list(), eqFun.bound=NULL,
+           ineqFun=list(), ineqFun.lower=NULL, ineqFun.upper=NULL,
+           control=list())
Error in loadNamespace(name) : there is no package called ‘Rdonlp2’
> donlp2NLPControl(
+     iterma=4000, nstep=20, fnscale=1, report=FALSE, rep.freq=1,
+     tau0=1, tau=0.1, del0=1, epsx=1e-05, delmin=0.1 * del0,
+     epsdif=1e-08, nreset.multiplier=1, difftype=3, epsfcn=1e-16,
+     taubnd=1, hessian=FALSE, te0=TRUE, te1=FALSE, te2=FALSE,
+     te3=FALSE, silent=TRUE, intakt=TRUE)
Error in donlp2NLPControl(iterma = 4000, nstep = 20, fnscale = 1, report = FALSE,  :
> rdonlp2
function (...)
{
Rdonlp2::donlp2(...)
}
<bytecode: 0x000000001b0c9470>
<environment: namespace:fPortfolio>
>
> rsolnpNLP(start, objective,
+           lower=0, upper=1, linCons, funCons, control=list())
Error in rsolnpNLP(start, objective, lower = 0, upper = 1, linCons, funCons,  :
> solnpNLP(start, objective,
+          par.lower=NULL, par.upper=NULL,
+          eqA=NULL, eqA.bound=NULL,
+          ineqA=NULL, ineqA.lower=NULL, ineqA.upper=NULL,
+          eqFun=list(), eqFun.bound=NULL,
+          ineqFun=list(), ineqFun.lower=NULL, ineqFun.upper=NULL,
+          control=list())
Error in solnpNLP(start, objective, par.lower = NULL, par.upper = NULL,  :
> solnpNLPControl(
+     rho=1, outer.iter=400, inner.iter=800, delta=1e-07, tol=1e-08, trace=0)
\$`rho`
 1

\$outer.iter
 400

\$inner.iter
 800

\$delta
 1e-07

\$tol
 1e-08

\$trace
 0

> rsolnp
>
> rnlminb2NLP(start, objective,
+             lower=0, upper=1, linCons, funCons, control=list())
Error in rnlminb2NLP(start, objective, lower = 0, upper = 1, linCons,  :
> nlminb2NLP(start, objective,
+            par.lower=NULL, par.upper=NULL,
+            eqA=NULL, eqA.bound=NULL,
+            ineqA=NULL, ineqA.lower=NULL, ineqA.upper=NULL,
+            eqFun=list(), eqFun.bound=NULL,
+            ineqFun=list(), ineqFun.lower=NULL, ineqFun.upper=NULL,
+            control=list())
Error in loadNamespace(name) : there is no package called ‘Rnlminb2’
> nlminb2NLPControl(
+     eval.max=500, iter.max=400, trace=0, abs.tol=1e-20, rel.tol=1e-10,
+     x.tol=1.5e-08, step.min=2.2e-14, scale=1, R=1, beta.tol=1e-20)
\$`eval.max`
 500

\$iter.max
 400

\$trace
 0

\$abs.tol
 1e-20

\$rel.tol
 1e-10

\$x.tol
 1.5e-08

\$step.min
 2.2e-14

\$scale
 1

\$R
 1

\$beta.tol
 1e-20

> rnlminb2
function (...)
{
Rnlminb2::nlminb2(...)
}
<bytecode: 0x000000001d797450>
<environment: namespace:fPortfolio>
>
> ramplNLP(start, objective,
+          lower=0, upper=1, amplCons, control=list(), ...)
Error: '...' used in an incorrect context
> amplNLP()
 NA
> amplNLPControl(
+     solver="minos", project="ampl", trace=FALSE)
\$`solver`
 "minos"

\$project
 "ampl"

\$trace
 FALSE

> rdonlp2NLP(start, objective,
+            lower=0, upper=1, linCons, funCons, control=list())
Error in rdonlp2NLP(start, objective, lower = 0, upper = 1, linCons, funCons,  :
> donlp2NLP(start, objective,
+           par.lower=NULL, par.upper=NULL,
+           eqA=NULL, eqA.bound=NULL,
+           ineqA=NULL, ineqA.lower=NULL, ineqA.upper=NULL,
+           eqFun=list(), eqFun.bound=NULL,
+           ineqFun=list(), ineqFun.lower=NULL, ineqFun.upper=NULL,
+           control=list())
Error in loadNamespace(name) : there is no package called ‘Rdonlp2’
> donlp2NLPControl(
+     iterma=4000, nstep=20, fnscale=1, report=FALSE, rep.freq=1,
+     tau0=1, tau=0.1, del0=1, epsx=1e-05, delmin=0.1 * del0,
+     epsdif=1e-08, nreset.multiplier=1, difftype=3, epsfcn=1e-16,
+     taubnd=1, hessian=FALSE, te0=TRUE, te1=FALSE, te2=FALSE,
+     te3=FALSE, silent=TRUE, intakt=TRUE)
Error in donlp2NLPControl(iterma = 4000, nstep = 20, fnscale = 1, report = FALSE,  :
> rdonlp2
function (...)
{
Rdonlp2::donlp2(...)
}
<bytecode: 0x000000001b0c9470>
<environment: namespace:fPortfolio>
> nlminb2NLPControl(
+     eval.max=500, iter.max=400, trace=0, abs.tol=1e-20, rel.tol=1e-10,
+     x.tol=1.5e-08, step.min=2.2e-14, scale=1, R=1, beta.tol=1e-20)
\$`eval.max`
 500

\$iter.max
 400

\$trace
 0

\$abs.tol
 1e-20

\$rel.tol
 1e-10

\$x.tol
 1.5e-08

\$step.min
 2.2e-14

\$scale
 1

\$R
 1

\$beta.tol
 1e-20

> rnlminb2
function (...)
{
Rnlminb2::nlminb2(...)
}
<bytecode: 0x000000001d797450>
<environment: namespace:fPortfolio>
> library("tseries", lib.loc="~/R/win-library/3.5")

‘tseries’ version: 0.10-46

‘tseries’ is a package for time series
analysis and computational finance.

See ‘library(help="tseries")’ for
details.

Warning message:
package ‘tseries’ was built under R version 3.5.2
> library(tseries)
> stabilityAnalytics(index, method=c("turns", "drawdowns", "garch",
+                                    "riskmetrics", "bcp", "pcout"), ...)
Error: '...' used in an incorrect context
> backtestPlot(object, which="all", labels=TRUE, legend=TRUE, at=NULL,
+              format=NULL, cex=0.6, font=1, family="mono")
Error in backtestAssetsPlot(object, labels, legend, at, format) :
>
> backtestAssetsPlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestAssetsPlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
> backtestWeightsPlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestWeightsPlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
> backtestRebalancePlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestRebalancePlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
> backtestPortfolioPlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestPortfolioPlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
> backtestDrawdownPlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestDrawdownPlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
> backtestReportPlot(object, cex=0.6, font=1, family="mono")
> backtestStats(100, FUN = "rollingSigma", ...)
Error: '...' used in an incorrect context
> library("PortfolioAnalytics", lib.loc="~/R/win-library/3.5")

Attaching package: ‘zoo’

The following object is masked from ‘package:timeSeries’:

time<-

The following objects are masked from ‘package:base’:

as.Date, as.Date.numeric

Copyright (c) 2004-2018 Peter Carl and Brian G. Peterson, GPL-2 | GPL-3
https://github.com/braverock/PerformanceAnalytics

Attaching package: ‘PerformanceAnalytics’

The following objects are masked from ‘package:timeDate’:

kurtosis, skewness

The following object is masked from ‘package:graphics’:

legend

Warning messages:
1: package ‘PortfolioAnalytics’ was built under R version 3.5.2
2: package ‘zoo’ was built under R version 3.5.1
3: package ‘xts’ was built under R version 3.5.1
4: package ‘foreach’ was built under R version 3.5.2
5: package ‘PerformanceAnalytics’ was built under R version 3.5.2
> library(PortfolioAnalytics)
> data(indexes)
US Bonds US Equities Int'l Equities
1980-01-31  -0.0272      0.0610         0.0462
1980-02-29  -0.0669      0.0031        -0.0040
1980-03-31   0.0053     -0.0987        -0.1188
1980-04-30   0.0992      0.0429         0.0864
1980-05-31   0.0000      0.0562         0.0446
1980-06-30   0.0605      0.0296         0.0600
Commodities US Tbill Inflation
1980-01-31      0.0568   0.0104    0.0149
1980-02-29     -0.0093   0.0106    0.0146
1980-03-31     -0.1625   0.0121    0.0120
1980-04-30      0.0357   0.0137    0.0095
1980-05-31      0.0573   0.0106    0.0095
1980-06-30      0.0533   0.0066    0.0000
Warning message:
index class is Date, which does not support timezones.
Expected 'UTC' timezone, but indexTZ is ''
> LPP <- LPP2005REC[1:12, 1:4]
> colnames(LPP) <- abbreviate(colnames(LPP), 2)
> finCenter(LPP) <- "GMT"
> plot(LPP[, 1], type = "o", col = "steelblue",
+      main = "LPP", xlab = "2005", ylab = "Return")
> plot(LPP[, 1], at="auto", type = "o", col = "steelblue",
+      main = "LPP", xlab = "2005", ylab = "Return")
> plot(LPP[, 1:2], type = "o", col = "steelblue",
+      main = "LPP", xlab = "2005", ylab = "Return"
+ )
> plot(LPP[, 1], LPP[, 2], type = "p", col = "steelblue",
+      main = "LPP", xlab = "Return 1", ylab = "Return 2")
> LPP <- as.timeSeries(data(LPP2005REC))
> ZRH <- as.timeSeries(LPP[,"SPI"], zone = "Zurich", FinCenter = "Zurich")
> NYC <- as.timeSeries(LPP[,"LMI"], zone = "NewYork", FinCenter = "NewYork")
> finCenter(ZRH)
 "Zurich"
> finCenter(NYC)
 "NewYork"
> plot(ZRH, at="auto", type = "p", pch = 19, col = "blue")
> points(NYC, pch = 19, col = "red")
> finCenter(ZRH) <- "Zurich"
> finCenter(NYC) <- "NewYork"
> at <- unique(round(time(ZRH)))
> plot(ZRH, type = "p", pch = 19, col = "blue", format = "%b %d", at = at,
+      FinCenter = "GMT", xlab = "GMT", main = "ZRH - GMT")
> points(NYC, FinCenter = "GMT", pch = 19, col = "red")
> data(edhec)
> ret <- edhec[, 1:4]
> pspec <- portfolio.spec(assets=colnames(ret))
> pspec <- add.constraint(pspec, type="box", min=0.05, max=0.45)
+                         type="box",
+                         min=c(0.05, 0.10, 0.08, 0.06),
+                         max=c(0.45, 0.55, 0.35, 0.65))
> centroid.buckets {PortfolioAnalytics}
Error: unexpected '{' in "centroid.buckets {"
> centroid.buckets(buckets, simulations = 1000)
Error in centroid.buckets(buckets, simulations = 1000) :
> chart.EfficientFrontier(object, ...,
+                         match.col = "ES", n.portfolios = 25, xlim = NULL, ylim = NULL,
+                         cex.axis = 0.8, element.color = "darkgray", main = "Efficient Frontier",
+                         RAR.text = "SR", rf = 0, tangent.line = TRUE, cex.legend = 0.8,
+                         chart.assets = TRUE, labels.assets = TRUE, pch.assets = 21,
+                         cex.assets = 0.8)
Error: '...' used in an incorrect context
> var.portfolio(R, weights)
Error in as.vector(x, mode) :
cannot coerce type 'closure' to vector of type 'any'
> Data
 NA
> library("FRAPO", lib.loc="~/R/win-library/3.5")
Using the GLPK callable library version 4.47
Financial Risk Modelling and Portfolio Optimisation with R (version 0.4-1)

Warning messages:
1: package ‘FRAPO’ was built under R version 3.5.2
2: package ‘cccp’ was built under R version 3.5.2
3: package ‘Rglpk’ was built under R version 3.5.2
4: package ‘slam’ was built under R version 3.5.2
> library(FRAPO)
> data(StockIndex)
> y <- StockIndex[, "SP500"]
> cs <- capser(y, min = 100, max = 200)
 200 200 200 200 200 200
> yret <- diff(log(y))
> bilson <- trdbilson(yret, exponent = 2)
  0.01948730 -0.01935785  0.01177193 -0.04517835
  0.10349955 -0.01545405
> es <- trdes(yret, lambda = 0.95)
  0.018485339 -0.017438597  0.010304526
 -0.042136563  0.093941299 -0.009969699
> cs <- capser(y, min = 100, max = 200)
> yret <- diff(log(y))
> binary <- trdbinary(yret)
  0.02477189 -0.02460780  0.01497858 -0.05712579
  0.12829276 -0.01965560
> hp <- trdhp(y, lambda = 1600)
 387.1268 389.8778 392.6293 395.3850 398.1460
 400.9114
> sma <- trdsma(y, n.periods = 24)
       NA       NA       NA       NA       NA
       NA       NA       NA       NA       NA
       NA       NA       NA       NA       NA
       NA       NA       NA       NA       NA
       NA       NA       NA 418.2483 420.7617
 423.6004 426.5617 429.7025 433.3096 435.5387
> wma <- trdwma(y, weights = c(0.4, 0.3, 0.2, 0.1))
      NA      NA      NA 391.205 384.938
 395.898 402.343 408.264 408.556 410.505
 412.633 411.220 416.691 416.362 417.228
 418.201 423.188 430.449 435.296 439.885
 445.320 444.131 446.805 448.474 448.468
 454.988 457.319 462.336 463.207 465.452
> require(graphics)
> x <- c(0,0,0,100,0,0,0)
> y <- c(0,0,1, 2 ,1,0,0)/4
> zapsmall(convolve(x, y))
 50 25  0  0  0  0 25
> zapsmall(convolve(x, y[3:5], type = "f"))
  0 25 50 25  0
> x <- rnorm(50)
> y <- rnorm(50)
> all.equal(convolve(x, y, conj = FALSE), rev(convolve(rev(y),x)))
 TRUE
> n <- length(x <- -20:24)
> y <- (x-10)^2/1000 + rnorm(x)/8
> Han <- function(y) # Hanning
+     convolve(y, c(1,2,1)/4, type = "filter")
> plot(x, y, main = "Using  convolve(.) for Hanning filters")
> lines(x[-c(1  , n)      ], Han(y), col = "red")
> lines(x[-c(1:2, (n-1):n)], Han(Han(y)), lwd = 2, col = "dark blue")
> library("cvar", lib.loc="~/R/win-library/3.5")

Attaching package: ‘cvar’

The following objects are masked from ‘package:PerformanceAnalytics’:

ES, VaR

Warning message:
package ‘cvar’ was built under R version 3.5.2
> library(cvar)
> mo <- GarchModel(omega = 0.4, alpha = 0.3, beta = 0.5)
> a1 <- sim_garch1c1(mo, n = 1000, n.start = 100)
> a.pred <- predict(mo, n.ahead = 5, Nsim = 1000, eps = a1\$eps, sigmasq = a1\$h, seed = 1234)
> a.pred\$eps
 0 0 0 0 0
> a.pred\$pi_plugin
lwr      upr
[1,] -1.992935 1.992935
[2,] -2.171179 2.171179
[3,] -2.303866 2.303866
[4,] -2.404750 2.404750
[5,] -2.482507 2.482507
> a.pred\$pi_sim
2.5%    97.5%
[1,] -1.962977 1.978131
[2,] -2.255450 2.206946
[3,] -2.182829 2.342846
[4,] -2.310415 2.150915
[5,] -2.440387 2.461284
> t(apply(a.pred\$dist_sim\$eps, 1, function(x) quantile(x, c(0.025, 0.975))))
2.5%    97.5%
[1,] -1.962977 1.978131
[2,] -2.255450 2.206946
[3,] -2.182829 2.342846
[4,] -2.310415 2.150915
[5,] -2.440387 2.461284
> t(apply(a.pred\$dist_sim\$eps, 1, function(x) summary(x)))
Min.    1st Qu.       Median
[1,] -3.146715 -0.6712101 -0.051807252
[2,] -5.517102 -0.7102704  0.046509597
[3,] -4.474419 -0.7719921 -0.013547190
[4,] -3.533618 -0.7864974 -0.001579576
[5,] -4.618284 -0.7965248 -0.018277768
Mean   3rd Qu.     Max.
[1,] -0.003854556 0.6966841 3.249664
[2,]  0.012472968 0.7621015 3.638396
[3,] -0.005678992 0.7372569 3.750280
[4,] -0.027251506 0.7057924 3.989376
[5,] -0.015917251 0.7528478 4.066457
> plot(density(a.pred\$dist_sim\$eps[2, ]), ylim = c(0,.25))
> lines(density(a.pred\$dist_sim\$eps[5, ]), col = "blue")
> a.pred\$h
 1.033928 1.227142 1.381714 1.505371 1.604297
> sqrt(a.pred\$h)
 1.016823 1.107765 1.175463 1.226936 1.266608
> plot(density(sqrt(a.pred\$dist_sim\$h[2, ])), xlim = c(0, 6))
> x <- seq(0.01, 0.99, length = 100)
> y <- sapply(x, function(p) cvar::ES(qnorm, x = p, dist.type = "qf"))
> yS <- sapply(x, function(p) - VaRES::esnormal(p))
Error in loadNamespace(name) : there is no package called ‘VaRES’
Called from: value[[3L]](cond)
Browse> plot(x, y)
Browse>
Found more than one class "xts" in cache; using the first, from namespace 'quantmod'
Also defined by ‘FRAPO’'
Also defined by ‘FRAPO’
> lines(x, yS, col = "blue")
> mo1a <- GarchModel(omega = 1, alpha = 0.3, beta = 0.5)
> mo1b <- GarchModel(omega = 1, alpha = 0.3, beta = 0.5, cond.dist = "norm")
>
>
> mo2 <- GarchModel(omega = 1, alpha = 0.3, beta = 0.5, esp0 = - 1.5, h0 = 4.96)
> mot <- GarchModel(omega = 1, alpha = 0.3, beta = 0.5, cond.dist = list("std", nu = 5)
+ )
> View(mot)
>
> View(mo1b)
> View(mo2)
> mo2
\$`omega`
 1

\$alpha
 0.3

\$beta
 0.5

\$esp0
 -1.5

\$h0
 4.96

attr(,"class")
 "garch1c1"
> mo <- GarchModel(omega = 0.4, alpha = 0.3, beta = 0.5)
> a1 <- sim_garch1c1(mo, n = 1000, n.start = 100)
> a.pred <- predict(mo, n.ahead = 5, Nsim = 1000, eps = a1\$eps, sigmasq = a1\$h, seed = 1234)
> a.pred\$eps
 0 0 0 0 0
> a.pred\$pi_plugin
lwr      upr
[1,] -3.438538 3.438538
[2,] -3.315934 3.315934
[3,] -3.214486 3.214486
[4,] -3.130961 3.130961
[5,] -3.062502 3.062502
> a.pred\$pi_sim
2.5%    97.5%
[1,] -3.386849 3.412995
[2,] -3.497097 3.303965
[3,] -2.941401 3.237982
[4,] -3.082723 2.933992
[5,] -3.021358 3.018965
> t(apply(a.pred\$dist_sim\$eps, 1, function(x) quantile(x, c(0.025, 0.975))
+ )
+ )
2.5%    97.5%
[1,] -3.386849 3.412995
[2,] -3.497097 3.303965
[3,] -2.941401 3.237982
[4,] -3.082723 2.933992
[5,] -3.021358 3.018965
> t(apply(a.pred\$dist_sim\$eps, 1, function(x) summary(x)))
Min.    1st Qu.       Median
[1,] -5.429227 -1.1580814 -0.089386336
[2,] -9.204317 -1.0739463  0.071430081
[3,] -6.710671 -1.0521054 -0.019127552
[4,] -4.758922 -0.9992271 -0.001989747
[5,] -6.334744 -0.9571489 -0.020867447
Mean   3rd Qu.     Max.
[1,] -0.006650510 1.2020332 5.606852
[2,]  0.014196985 1.1367528 5.709787
[3,] -0.002970354 0.9980880 5.437128
[4,] -0.034739711 0.8886597 4.829568
[5,] -0.017372954 0.8908675 6.468912
> plot(density(a.pred\$dist_sim\$eps[2, ]), ylim = c(0,.25))
> lines(density(a.pred\$dist_sim\$eps[5, ]), col = "blue")
> cvar::VaR(qnorm, x = c(0.01, 0.05), dist.type = "qf")
 2.326348 1.644854
> muA <- 0.006408553
> sigma2A <- 0.0004018977
> res1 <- cvar::VaR(qnorm, x = 0.05, mean = muA, sd = sqrt(sigma2A))
> res2 <- cvar::VaR(qnorm, x = 0.05, intercept = muA, slope = sqrt(sigma2A))
> abs((res2 - res1)) # 0, intercept/slope equivalent to mean/sd
 0
> res1a <- cvar::VaR(pnorm, x = 0.05, dist.type = "cdf", mean = muA, sd = sqrt(sigma2A))
> res2a <- cvar::VaR(pnorm, x = 0.05, dist.type = "cdf", intercept = muA, slope = sqrt(sigma2A))
> abs((res1a - res2)) # 3.287939e-09
 3.287939e-09
> abs((res2a - res2)) # 5.331195e-11, intercept/slope better numerically
 5.331195e-11
> res1b <- cvar::VaR(pnorm, x = 0.05, dist.type = "cdf",
+                    mean = muA, sd = sqrt(sigma2A), tol = .Machine\$double.eps^0.75)
> res2b <- cvar::VaR(pnorm, x = 0.05, dist.type = "cdf",
+                    intercept = muA, slope = sqrt(sigma2A), tol = .Machine\$double.eps^0.75)
> abs((res1b - res2)) # 6.938894e-18 # both within machine precision
 6.938894e-18
> abs((res2b - res2)) # 1.040834e-16
 1.040834e-16
> abs((res1b - res2)/res2) # 2.6119e-16 # both within machine precision
 2.6119e-16
> abs((res2b - res2)/res2) # 3.91785e-15
 3.91785e-15
> if (requireNamespace("PerformanceAnalytics", quietly = TRUE)) withAutoprint({
+     data(edhec, package = "PerformanceAnalytics")
+     mu <- apply(edhec, 2, mean)
+     sigma2 <- apply(edhec, 2, var)
+     musigma2 <- cbind(mu, sigma2)
+ ## analogous calc. with PerformanceAnalytics::VaR
+ vPA <- apply(musigma2, 1, function(x)
+     PerformanceAnalytics::VaR(p = .95, method = "gaussian", invert = FALSE,
+                               mu = x, sigma = x, weights = 1)
+ max(abs((vPA - vAz1))) # 5.551115e-17
Error: unexpected symbol in:
"                              mu = x, sigma = x, weights = 1)
max"
> max(abs((vPA - vAz2))) #   ""
>
> max(abs((vPA - vAz1a))) # 3.287941e-09
> max(abs((vPA - vAz2a))) #  1.465251e-10, intercept/slope better
>
> max(abs((vPA - vAz1b))) # 4.374869e-13
> max(abs((vPA - vAz2b))) # 3.330669e-16
> }
Error: unexpected '}' in "}"
> library("scenario", lib.loc="~/R/win-library/3.5")
Warning message:
package ‘scenario’ was built under R version 3.5.2
> library(scenario)
> centroids <- cbind(c(0,2,3), c(0,2,1), c(0,-2,-3),c(0,-2,-1))
There were 36 warnings (use warnings() to see them)
> matplot(centroids, type="l", lwd = 3, col = "black", lty = 3)
> scenarios <- matrix(rep(centroids,5), ncol=20) + matrix(rnorm(60,0,0.25),ncol=20)
> matlines(scenarios, col = "grey")
> treeStruct <- rbind(c(1,1,1,1),
+                     c(2,2,5,5),
+                     c(3,4,6,7))
> checktree(treeStruct)
>

> tree <- buildtree(scenarios, treeStruct, jMax = 1000)
> matlines(centroids,lwd = 3, col = "black", lty = 3)
> matlines(centroids,lwd = 3, col = "green", lty = 3)

## Sunday, 27 January 2019

### R First attempt at a vector field

legend.position=c(1, 0.55), # Put legend inside plot area
+           legend.justification=c(1, 0.5))
> ggplot(tophit, aes(x=avg, y=name)) +
+     geom_segment(aes(yend=name), xend=0, colour="grey50") +
+     geom_point(size=3, aes(colour=lg)) +
+     scale_colour_brewer(palette="Set1", limits=c("NL","AL"), guide=FALSE) +
+     theme_bw() +
+     theme(panel.grid.major.y = element_blank()) +
+     facet_grid(lg ~ ., scales="free_y", space="free_y")
>
> ggplot(worldpop, aes(x=Year, y=Population)) + geom_line() + geom_point()
>
> ggplot(worldpop, aes(x=Year, y=Population)) + geom_line() + geom_point() +
+     scale_y_log10()
>
> tg <- ddply(ToothGrowth, c("supp", "dose"), summarise, length=mean(len))
> ggplot(tg, aes(x=dose, y=length, colour=supp)) + geom_line()
> ggplot(tg, aes(x=dose, y=length, linetype=supp)) + geom_line()
> sunspotyear <- data.frame(
+     Year = as.numeric(time(sunspot.year)),
+     Sunspots = as.numeric(sunspot.year)
+ )
>
> ggplot(sunspotyear, aes(x=Year, y=Sunspots)) + geom_area()
> ggplot(sunspotyear, aes(x=Year, y=Sunspots)) +
+     geom_area(colour="black", fill="blue", alpha=.2)
> ggplot(sunspotyear, aes(x=Year, y=Sunspots)) +
+     geom_area(fill="blue", alpha=.2) +
+     geom_line()
>
> ggplot(uspopage, aes(x=Year, y=Thousands, fill=AgeGroup)) + geom_area()
> ggplot(uspopage, aes(x=Year, y=Thousands, fill=AgeGroup)) +
+     geom_area(colour="black", size=.2, alpha=.4) +
+     scale_fill_brewer(palette="Blues", breaks=rev(levels(uspopage\$AgeGroup)))
>
> ggplot(uspopage, aes(x=Year, y=Thousands, fill=AgeGroup, order=desc(AgeGroup))) +
+     geom_area(colour=NA, alpha=.4) +
+     scale_fill_brewer(palette="Blues") +
+     geom_line(position="stack", size=.2)
>
> uspopage_prop <- ddply(uspopage, "Year", transform,
+                        Percent = Thousands / sum(Thousands) * 100)
>
> ggplot(uspopage_prop, aes(x=Year, y=Percent, fill=AgeGroup)) +
+     geom_area(colour="black", size=.2, alpha=.4) +
+     scale_fill_brewer(palette="Blues", breaks=rev(levels(uspopage\$AgeGroup)))
>
> clim <- subset(climate, Source == "Berkeley",
+                select=c("Year", "Anomaly10y", "Unc10y"))
>
> ggplot(clim, aes(x=Year, y=Anomaly10y)) +
+     geom_ribbon(aes(ymin=Anomaly10y-Unc10y, ymax=Anomaly10y+Unc10y),
+                 alpha=0.2) +
+     geom_line()
>
> ggplot(clim, aes(x=Year, y=Anomaly10y)) +
+     geom_line(aes(y=Anomaly10y-Unc10y), colour="grey50", linetype="dotted") +
+     geom_line(aes(y=Anomaly10y+Unc10y), colour="grey50", linetype="dotted") +
+     geom_line()
>
> heightweight[, c("ageYear", "heightIn")]
ageYear heightIn
1     11.92     56.3
2     12.92     62.3
3     12.75     63.3
4     13.42     59.0
5     15.92     62.5
6     14.25     62.5
7     15.42     59.0
8     11.83     56.5
9     13.33     62.0
10    11.67     53.8
11    11.58     61.5
12    14.83     61.5
13    13.08     64.5
14    12.42     58.3
15    11.92     51.3
16    12.08     58.8
17    15.92     65.3
18    12.50     59.5
19    12.25     61.3
20    15.00     63.3
21    11.75     61.8
22    11.67     53.5
23    13.67     58.0
24    14.67     61.3
25    15.42     63.3
26    13.83     61.5
27    14.58     60.8
28    15.00     59.0
29    17.50     65.5
30    12.17     56.3
31    14.17     64.3
32    13.50     58.0
33    12.42     64.3
34    11.58     57.5
35    15.50     57.8
36    16.42     61.5
37    14.08     62.3
38    14.75     61.8
39    15.42     65.3
40    15.17     58.3
41    14.42     62.8
42    13.83     59.3
43    14.00     61.5
44    14.08     62.0
45    12.50     61.3
46    15.33     62.3
47    11.58     52.8
48    12.25     59.8
49    12.00     59.5
50    14.75     61.3
51    14.83     63.5
52    16.42     64.8
53    12.17     60.0
54    12.08     59.0
55    12.25     55.8
56    12.08     57.8
57    12.92     61.3
58    13.92     62.3
59    15.25     64.3
60    11.92     55.5
61    15.25     64.5
62    15.42     60.0
63    12.33     56.3
64    12.25     58.3
65    12.83     60.0
66    13.00     54.5
67    12.00     55.8
68    12.83     62.8
69    12.67     60.5
70    15.92     63.3
71    15.83     66.8
72    11.67     60.0
73    12.33     60.5
74    15.75     64.3
75    11.92     58.3
76    14.83     66.5
77    13.67     65.3
78    13.08     60.5
79    12.25     59.5
80    12.33     59.0
81    14.75     61.3
82    14.25     61.5
83    14.33     64.8
84    15.83     56.8
85    15.25     66.5
86    11.92     61.5
87    14.92     63.0
88    15.50     57.0
89    15.17     65.5
90    15.17     62.0
91    11.83     56.0
92    13.75     61.3
93    13.75     55.5
94    12.83     61.0
95    12.50     54.5
96    12.92     66.0
97    13.58     56.5
98    11.75     56.0
99    12.25     51.5
100   17.50     62.0
101   14.25     63.0
102   13.92     61.0
103   15.17     64.0
104   12.00     61.0
105   16.08     59.8
106   11.75     61.3
107   13.67     63.3
108   15.50     63.5
109   14.08     61.5
110   14.58     60.3
111   15.00     61.3
112   13.75     64.8
113   13.08     60.5
114   12.00     57.3
115   12.50     59.5
116   12.50     60.8
117   11.58     60.5
118   15.75     67.0
119   15.25     64.8
120   12.25     50.5
121   12.17     57.5
122   13.33     60.5
123   13.00     61.8
124   14.42     61.3
125   12.58     66.3
126   11.75     53.3
127   12.50     59.0
128   13.67     57.8
129   12.75     60.0
130   17.17     68.3
132   14.67     63.8
133   14.67     65.0
134   11.67     59.5
135   15.42     66.0
136   15.00     61.8
137   12.17     57.3
138   15.25     66.0
139   11.67     56.5
140   12.58     58.3
141   12.58     61.0
142   12.00     62.8
143   13.33     59.3
144   14.83     67.3
145   16.08     66.3
146   13.50     64.5
147   13.67     60.5
148   15.50     66.0
149   11.92     57.5
150   14.58     64.0
151   14.58     68.0
152   14.58     63.5
153   14.42     69.0
154   14.17     63.8
155   14.50     66.0
156   13.67     63.5
157   12.00     59.5
158   13.00     66.3
159   12.42     57.0
160   12.00     60.0
161   12.25     57.0
162   15.67     67.3
163   14.08     62.0
164   14.33     65.0
165   12.50     59.5
166   16.08     67.8
167   13.08     58.0
168   14.00     60.0
169   11.67     58.5
170   13.00     58.3
171   13.00     61.5
172   13.17     65.0
173   15.33     66.5
174   13.00     68.5
175   12.00     57.0
176   14.67     61.5
177   14.00     66.5
178   12.42     52.5
179   11.83     55.0
180   15.67     71.0
181   16.92     66.5
182   11.83     58.8
183   15.75     66.3
184   15.67     65.8
185   16.67     71.0
186   12.67     59.5
187   14.50     69.8
188   13.83     62.5
189   12.08     56.5
190   11.92     57.5
191   13.58     65.3
192   13.83     67.3
193   15.17     67.0
194   14.42     66.0
195   12.92     61.8
196   13.50     60.0
197   14.75     63.0
198   14.75     60.5
199   14.58     65.5
200   13.83     62.0
201   12.50     59.0
202   12.50     61.8
203   15.67     63.3
204   13.58     66.0
205   14.25     61.8
206   13.50     63.0
207   11.75     57.5
208   14.50     63.0
209   11.83     56.0
210   12.33     60.5
211   11.67     56.8
212   13.33     64.0
213   12.00     60.0
214   17.17     69.5
215   13.25     63.3
216   12.42     56.3
217   16.08     72.0
218   16.17     65.3
219   12.67     60.8
220   12.17     55.0
221   11.58     55.0
222   15.50     66.5
223   13.42     56.8
224   12.75     64.8
225   16.33     64.5
226   13.67     58.0
227   13.25     62.8
228   14.83     63.8
229   12.75     57.8
230   12.92     57.3
231   14.83     63.5
232   11.83     55.0
233   13.67     66.5
234   15.75     65.0
235   13.67     61.5
236   13.92     62.0
237   12.58     59.3
>
> ggplot(heightweight, aes(x=ageYear, y=heightIn)) + geom_point()
>
> ggplot(heightweight, aes(x=ageYear, y=heightIn)) + geom_point(shape=21)
>
> ggplot(heightweight, aes(x=ageYear, y=heightIn)) + geom_point(size=1.5)
> ggplot(heightweight, aes(x=ageYear, y=heightIn)) + geom_point(size=1.5)
> heightweight[, c("sex", "ageYear", "heightIn")]
sex ageYear heightIn
1     f   11.92     56.3
2     f   12.92     62.3
3     f   12.75     63.3
4     f   13.42     59.0
5     f   15.92     62.5
6     f   14.25     62.5
7     f   15.42     59.0
8     f   11.83     56.5
9     f   13.33     62.0
10    f   11.67     53.8
11    f   11.58     61.5
12    f   14.83     61.5
13    f   13.08     64.5
14    f   12.42     58.3
15    f   11.92     51.3
16    f   12.08     58.8
17    f   15.92     65.3
18    f   12.50     59.5
19    f   12.25     61.3
20    f   15.00     63.3
21    f   11.75     61.8
22    f   11.67     53.5
23    f   13.67     58.0
24    f   14.67     61.3
25    f   15.42     63.3
26    f   13.83     61.5
27    f   14.58     60.8
28    f   15.00     59.0
29    f   17.50     65.5
30    f   12.17     56.3
31    f   14.17     64.3
32    f   13.50     58.0
33    f   12.42     64.3
34    f   11.58     57.5
35    f   15.50     57.8
36    f   16.42     61.5
37    f   14.08     62.3
38    f   14.75     61.8
39    f   15.42     65.3
40    f   15.17     58.3
41    f   14.42     62.8
42    f   13.83     59.3
43    f   14.00     61.5
44    f   14.08     62.0
45    f   12.50     61.3
46    f   15.33     62.3
47    f   11.58     52.8
48    f   12.25     59.8
49    f   12.00     59.5
50    f   14.75     61.3
51    f   14.83     63.5
52    f   16.42     64.8
53    f   12.17     60.0
54    f   12.08     59.0
55    f   12.25     55.8
56    f   12.08     57.8
57    f   12.92     61.3
58    f   13.92     62.3
59    f   15.25     64.3
60    f   11.92     55.5
61    f   15.25     64.5
62    f   15.42     60.0
63    f   12.33     56.3
64    f   12.25     58.3
65    f   12.83     60.0
66    f   13.00     54.5
67    f   12.00     55.8
68    f   12.83     62.8
69    f   12.67     60.5
70    f   15.92     63.3
71    f   15.83     66.8
72    f   11.67     60.0
73    f   12.33     60.5
74    f   15.75     64.3
75    f   11.92     58.3
76    f   14.83     66.5
77    f   13.67     65.3
78    f   13.08     60.5
79    f   12.25     59.5
80    f   12.33     59.0
81    f   14.75     61.3
82    f   14.25     61.5
83    f   14.33     64.8
84    f   15.83     56.8
85    f   15.25     66.5
86    f   11.92     61.5
87    f   14.92     63.0
88    f   15.50     57.0
89    f   15.17     65.5
90    f   15.17     62.0
91    f   11.83     56.0
92    f   13.75     61.3
93    f   13.75     55.5
94    f   12.83     61.0
95    f   12.50     54.5
96    f   12.92     66.0
97    f   13.58     56.5
98    f   11.75     56.0
99    f   12.25     51.5
100   f   17.50     62.0
101   f   14.25     63.0
102   f   13.92     61.0
103   f   15.17     64.0
104   f   12.00     61.0
105   f   16.08     59.8
106   f   11.75     61.3
107   f   13.67     63.3
108   f   15.50     63.5
109   f   14.08     61.5
110   f   14.58     60.3
111   f   15.00     61.3
112   m   13.75     64.8
113   m   13.08     60.5
114   m   12.00     57.3
115   m   12.50     59.5
116   m   12.50     60.8
117   m   11.58     60.5
118   m   15.75     67.0
119   m   15.25     64.8
120   m   12.25     50.5
121   m   12.17     57.5
122   m   13.33     60.5
123   m   13.00     61.8
124   m   14.42     61.3
125   m   12.58     66.3
126   m   11.75     53.3
127   m   12.50     59.0
128   m   13.67     57.8
129   m   12.75     60.0
130   m   17.17     68.3
132   m   14.67     63.8
133   m   14.67     65.0
134   m   11.67     59.5
135   m   15.42     66.0
136   m   15.00     61.8
137   m   12.17     57.3
138   m   15.25     66.0
139   m   11.67     56.5
140   m   12.58     58.3
141   m   12.58     61.0
142   m   12.00     62.8
143   m   13.33     59.3
144   m   14.83     67.3
145   m   16.08     66.3
146   m   13.50     64.5
147   m   13.67     60.5
148   m   15.50     66.0
149   m   11.92     57.5
150   m   14.58     64.0
151   m   14.58     68.0
152   m   14.58     63.5
153   m   14.42     69.0
154   m   14.17     63.8
155   m   14.50     66.0
156   m   13.67     63.5
157   m   12.00     59.5
158   m   13.00     66.3
159   m   12.42     57.0
160   m   12.00     60.0
161   m   12.25     57.0
162   m   15.67     67.3
163   m   14.08     62.0
164   m   14.33     65.0
165   m   12.50     59.5
166   m   16.08     67.8
167   m   13.08     58.0
168   m   14.00     60.0
169   m   11.67     58.5
170   m   13.00     58.3
171   m   13.00     61.5
172   m   13.17     65.0
173   m   15.33     66.5
174   m   13.00     68.5
175   m   12.00     57.0
176   m   14.67     61.5
177   m   14.00     66.5
178   m   12.42     52.5
179   m   11.83     55.0
180   m   15.67     71.0
181   m   16.92     66.5
182   m   11.83     58.8
183   m   15.75     66.3
184   m   15.67     65.8
185   m   16.67     71.0
186   m   12.67     59.5
187   m   14.50     69.8
188   m   13.83     62.5
189   m   12.08     56.5
190   m   11.92     57.5
191   m   13.58     65.3
192   m   13.83     67.3
193   m   15.17     67.0
194   m   14.42     66.0
195   m   12.92     61.8
196   m   13.50     60.0
197   m   14.75     63.0
198   m   14.75     60.5
199   m   14.58     65.5
200   m   13.83     62.0
201   m   12.50     59.0
202   m   12.50     61.8
203   m   15.67     63.3
204   m   13.58     66.0
205   m   14.25     61.8
206   m   13.50     63.0
207   m   11.75     57.5
208   m   14.50     63.0
209   m   11.83     56.0
210   m   12.33     60.5
211   m   11.67     56.8
212   m   13.33     64.0
213   m   12.00     60.0
214   m   17.17     69.5
215   m   13.25     63.3
216   m   12.42     56.3
217   m   16.08     72.0
218   m   16.17     65.3
219   m   12.67     60.8
220   m   12.17     55.0
221   m   11.58     55.0
222   m   15.50     66.5
223   m   13.42     56.8
224   m   12.75     64.8
225   m   16.33     64.5
226   m   13.67     58.0
227   m   13.25     62.8
228   m   14.83     63.8
229   m   12.75     57.8
230   m   12.92     57.3
231   m   14.83     63.5
232   m   11.83     55.0
233   m   13.67     66.5
234   m   15.75     65.0
235   m   13.67     61.5
236   m   13.92     62.0
237   m   12.58     59.3
> ggplot(heightweight, aes(x=ageYear, y=heightIn, colour=sex)) + geom_point()
>
> ggplot(heightweight, aes(x=ageYear, y=heightIn, shape=sex)) + geom_point()
> ggplot(heightweight, aes(x=ageYear, y=heightIn, shape=sex, colour=sex)) +
+     geom_point()
> ggplot(heightweight, aes(x=ageYear, y=heightIn, shape=sex, colour=sex)) +
+     geom_point() +
+     scale_shape_manual(values=c(1,2)) +
+     scale_colour_brewer(palette="Set1")
> ggplot(heightweight, aes(x=ageYear, y=heightIn)) + geom_point(shape=3)
>
> ggplot(heightweight, aes(x=ageYear, y=heightIn, shape=sex)) +
+     geom_point(size=3) + scale_shape_manual(values=c(1, 4))
> hw\$weightGroup <- cut(hw\$weightLb, breaks=c(-Inf, 100, Inf),
+                       labels=c("< 100", ">= 100"))
Error in cut(hw\$weightLb, breaks = c(-Inf, 100, Inf), labels = c("< 100",  :
>
> hw <- heightweight
> hw\$weightGroup <- cut(hw\$weightLb, breaks=c(-Inf, 100, Inf),
+                       labels=c("< 100", ">= 100"))
>
> ggplot(hw, aes(x=ageYear, y=heightIn, shape=sex, fill=weightGroup)) +
+     geom_point(size=2.5) +
+     scale_shape_manual(values=c(21, 24)) +
+     scale_fill_manual(values=c(NA, "black"),
+                       guide=guide_legend(override.aes=list(shape=21)))
> heightweight[, c("sex", "ageYear", "heightIn", "weightLb")]
sex ageYear heightIn weightLb
1     f   11.92     56.3     85.0
2     f   12.92     62.3    105.0
3     f   12.75     63.3    108.0
4     f   13.42     59.0     92.0
5     f   15.92     62.5    112.5
6     f   14.25     62.5    112.0
7     f   15.42     59.0    104.0
8     f   11.83     56.5     69.0
9     f   13.33     62.0     94.5
10    f   11.67     53.8     68.5
11    f   11.58     61.5    104.0
12    f   14.83     61.5    103.5
13    f   13.08     64.5    123.5
14    f   12.42     58.3     93.0
15    f   11.92     51.3     50.5
16    f   12.08     58.8     89.0
17    f   15.92     65.3    107.0
18    f   12.50     59.5     78.5
19    f   12.25     61.3    115.0
20    f   15.00     63.3    114.0
21    f   11.75     61.8     85.0
22    f   11.67     53.5     81.0
23    f   13.67     58.0     83.5
24    f   14.67     61.3    112.0
25    f   15.42     63.3    101.0
26    f   13.83     61.5    103.5
27    f   14.58     60.8     93.5
28    f   15.00     59.0    112.0
29    f   17.50     65.5    140.0
30    f   12.17     56.3     83.5
31    f   14.17     64.3     90.0
32    f   13.50     58.0     84.0
33    f   12.42     64.3    110.5
34    f   11.58     57.5     96.0
35    f   15.50     57.8     95.0
36    f   16.42     61.5    121.0
37    f   14.08     62.3     99.5
38    f   14.75     61.8    142.5
39    f   15.42     65.3    118.0
40    f   15.17     58.3    104.5
41    f   14.42     62.8    102.5
42    f   13.83     59.3     89.5
43    f   14.00     61.5     95.0
44    f   14.08     62.0     98.5
45    f   12.50     61.3     94.0
46    f   15.33     62.3    108.0
47    f   11.58     52.8     63.5
48    f   12.25     59.8     84.5
49    f   12.00     59.5     93.5
50    f   14.75     61.3    112.0
51    f   14.83     63.5    148.5
52    f   16.42     64.8    112.0
53    f   12.17     60.0    109.0
54    f   12.08     59.0     91.5
55    f   12.25     55.8     75.0
56    f   12.08     57.8     84.0
57    f   12.92     61.3    107.0
58    f   13.92     62.3     92.5
59    f   15.25     64.3    109.5
60    f   11.92     55.5     84.0
61    f   15.25     64.5    102.5
62    f   15.42     60.0    106.0
63    f   12.33     56.3     77.0
64    f   12.25     58.3    111.5
65    f   12.83     60.0    114.0
66    f   13.00     54.5     75.0
67    f   12.00     55.8     73.5
68    f   12.83     62.8     93.5
69    f   12.67     60.5    105.0
70    f   15.92     63.3    113.5
71    f   15.83     66.8    140.0
72    f   11.67     60.0     77.0
73    f   12.33     60.5     84.5
74    f   15.75     64.3    113.5
75    f   11.92     58.3     77.5
76    f   14.83     66.5    117.5
77    f   13.67     65.3     98.0
78    f   13.08     60.5    112.0
79    f   12.25     59.5    101.0
80    f   12.33     59.0     95.0
81    f   14.75     61.3     81.0
82    f   14.25     61.5     91.0
83    f   14.33     64.8    142.0
84    f   15.83     56.8     98.5
85    f   15.25     66.5    112.0
86    f   11.92     61.5    116.5
87    f   14.92     63.0     98.5
88    f   15.50     57.0     83.5
89    f   15.17     65.5    133.0
90    f   15.17     62.0     91.5
91    f   11.83     56.0     72.5
92    f   13.75     61.3    106.5
93    f   13.75     55.5     67.0
94    f   12.83     61.0    122.5
95    f   12.50     54.5     74.0
96    f   12.92     66.0    144.5
97    f   13.58     56.5     84.0
98    f   11.75     56.0     72.5
99    f   12.25     51.5     64.0
100   f   17.50     62.0    116.0
101   f   14.25     63.0     84.0
102   f   13.92     61.0     93.5
103   f   15.17     64.0    111.5
104   f   12.00     61.0     92.0
105   f   16.08     59.8    115.0
106   f   11.75     61.3     85.0
107   f   13.67     63.3    108.0
108   f   15.50     63.5    108.0
109   f   14.08     61.5     85.0
110   f   14.58     60.3     86.0
111   f   15.00     61.3    110.5
112   m   13.75     64.8     98.0
113   m   13.08     60.5    105.0
114   m   12.00     57.3     76.5
115   m   12.50     59.5     84.0
116   m   12.50     60.8    128.0
117   m   11.58     60.5     87.0
118   m   15.75     67.0    128.0
119   m   15.25     64.8    111.0
120   m   12.25     50.5     79.0
121   m   12.17     57.5     90.0
122   m   13.33     60.5     84.0
123   m   13.00     61.8    112.0
124   m   14.42     61.3     93.0
125   m   12.58     66.3    117.0
126   m   11.75     53.3     84.0
127   m   12.50     59.0     99.5
128   m   13.67     57.8     95.0
129   m   12.75     60.0     84.0
130   m   17.17     68.3    134.0
132   m   14.67     63.8     98.5
133   m   14.67     65.0    118.5
134   m   11.67     59.5     94.5
135   m   15.42     66.0    105.0
136   m   15.00     61.8    104.0
137   m   12.17     57.3     83.0
138   m   15.25     66.0    105.5
139   m   11.67     56.5     84.0
140   m   12.58     58.3     86.0
141   m   12.58     61.0     81.0
142   m   12.00     62.8     94.0
143   m   13.33     59.3     78.5
144   m   14.83     67.3    119.5
145   m   16.08     66.3    133.0
146   m   13.50     64.5    119.0
147   m   13.67     60.5     95.0
148   m   15.50     66.0    112.0
149   m   11.92     57.5     75.0
150   m   14.58     64.0     92.0
151   m   14.58     68.0    112.0
152   m   14.58     63.5     98.5
153   m   14.42     69.0    112.5
154   m   14.17     63.8    112.5
155   m   14.50     66.0    108.0
156   m   13.67     63.5    108.0
157   m   12.00     59.5     88.0
158   m   13.00     66.3    106.0
159   m   12.42     57.0     92.0
160   m   12.00     60.0    117.5
161   m   12.25     57.0     84.0
162   m   15.67     67.3    112.0
163   m   14.08     62.0    100.0
164   m   14.33     65.0    112.0
165   m   12.50     59.5     84.0
166   m   16.08     67.8    127.5
167   m   13.08     58.0     80.5
168   m   14.00     60.0     93.5
169   m   11.67     58.5     86.5
170   m   13.00     58.3     92.5
171   m   13.00     61.5    108.5
172   m   13.17     65.0    121.0
173   m   15.33     66.5    112.0
174   m   13.00     68.5    114.0
175   m   12.00     57.0     84.0
176   m   14.67     61.5     81.0
177   m   14.00     66.5    111.5
178   m   12.42     52.5     81.0
179   m   11.83     55.0     70.0
180   m   15.67     71.0    140.0
181   m   16.92     66.5    117.0
182   m   11.83     58.8     84.0
183   m   15.75     66.3    112.0
184   m   15.67     65.8    150.5
185   m   16.67     71.0    147.0
186   m   12.67     59.5    105.0
187   m   14.50     69.8    119.5
188   m   13.83     62.5     84.0
189   m   12.08     56.5     91.0
190   m   11.92     57.5    101.0
191   m   13.58     65.3    117.5
192   m   13.83     67.3    121.0
193   m   15.17     67.0    133.0
194   m   14.42     66.0    112.0
195   m   12.92     61.8     91.5
196   m   13.50     60.0    105.0
197   m   14.75     63.0    111.0
198   m   14.75     60.5    112.0
199   m   14.58     65.5    114.0
200   m   13.83     62.0     91.0
201   m   12.50     59.0     98.0
202   m   12.50     61.8    118.0
203   m   15.67     63.3    115.5
204   m   13.58     66.0    112.0
205   m   14.25     61.8    112.0
206   m   13.50     63.0     91.0
207   m   11.75     57.5     85.0
208   m   14.50     63.0    112.0
209   m   11.83     56.0     87.5
210   m   12.33     60.5    118.0
211   m   11.67     56.8     83.5
212   m   13.33     64.0    116.0
213   m   12.00     60.0     89.0
214   m   17.17     69.5    171.5
215   m   13.25     63.3    112.0
216   m   12.42     56.3     72.0
217   m   16.08     72.0    150.0
218   m   16.17     65.3    134.5
219   m   12.67     60.8     97.0
220   m   12.17     55.0     71.5
221   m   11.58     55.0     73.5
222   m   15.50     66.5    112.0
223   m   13.42     56.8     75.0
224   m   12.75     64.8    128.0
225   m   16.33     64.5     98.0
226   m   13.67     58.0     84.0
227   m   13.25     62.8     99.0
228   m   14.83     63.8    112.0
229   m   12.75     57.8     79.5
230   m   12.92     57.3     80.5
231   m   14.83     63.5    102.5
232   m   11.83     55.0     76.0
233   m   13.67     66.5    112.0
234   m   15.75     65.0    114.0
235   m   13.67     61.5    140.0
236   m   13.92     62.0    107.5
237   m   12.58     59.3     87.0
>
> ggplot(heightweight, aes(x=ageYear, y=heightIn, colour=weightLb)) + geom_point()
> ggplot(heightweight, aes(x=ageYear, y=heightIn, size=weightLb)) + geom_point()
> ggplot(heightweight, aes(x=weightLb, y=heightIn, fill=ageYear)) +
+     geom_point(shape=21, size=2.5) +
>
> ggplot(heightweight, aes(x=weightLb, y=heightIn, fill=ageYear)) +
+     geom_point(shape=21, size=2.5) +
+                         guide=guide_legend())
>
> ggplot(heightweight, aes(x=ageYear, y=heightIn, size=weightLb, colour=sex)) +
+     geom_point(alpha=.5) +
+     scale_size_area() + # Make area proportional to numeric value
+     scale_colour_brewer(palette="Set1")
>
> sp <- ggplot(diamonds, aes(x=carat, y=price))
>
> sp + geom_point()
>
> sp + geom_point(alpha=.1)
> sp + geom_point(alpha=.01)
> sp + stat_bin2d()
>
> sp + stat_bin2d(bins=50) +
>
> library(hexbin)
Warning message:
package ‘hexbin’ was built under R version 3.5.2
> library(hexbin)
> sp + stat_binhex() +
+                         limits=c(0, 8000))
>
> sp + stat_binhex() +
+                         breaks=c(0, 250, 500, 1000, 2000, 4000, 6000),
+                         limits=c(0, 6000))
>
> sp1 <- ggplot(ChickWeight, aes(x=Time, y=weight))
>
> sp1 + geom_point()
>
> sp1 + geom_point(position="jitter")
> sp1 + geom_point(position=position_jitter(width=.5, height=0))
>
> sp1 + geom_boxplot(aes(group=Time))
> sp <- ggplot(heightweight, aes(x=ageYear, y=heightIn))
>
> sp + geom_point() + stat_smooth(method=lm)
> sp + geom_point() + stat_smooth(method=lm, level=0.99)
> sp + geom_point() + stat_smooth(method=lm, se=FALSE)
>
> sp + geom_point(colour="grey60") +
+     stat_smooth(method=lm, se=FALSE, colour="black")
>
> sp + geom_point(colour="grey60") + stat_smooth()
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
> sp + geom_point(colour="grey60") + stat_smooth(method=loess)
>
> library(MASS)
> b <- biopsy
> b\$classn[b\$class=="benign"] <- 0
> b\$classn[b\$class=="malignant"] <- 1
>
> ggplot(b, aes(x=V1, y=classn)) +
+     geom_point(position=position_jitter(width=0.3, height=0.06), alpha=0.4,
+                shape=21, size=1.5) +
+     stat_smooth(method=glm, family=binomial)

Warning: Ignoring unknown parameters: family
>
> sps <- ggplot(heightweight, aes(x=ageYear, y=heightIn, colour=sex)) +
+     geom_point() +
+     scale_colour_brewer(palette="Set1")
>
> sps + geom_smooth()
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
>
> sps + geom_smooth(method=lm, se=FALSE, fullrange=TRUE)
> model <- lm(heightIn ~ ageYear + I(ageYear^2), heightweight)
>
> xmin <- min(heightweight\$ageYear)
>
> xmax <- max(heightweight\$ageYear)
> predicted <- data.frame(ageYear=seq(xmin, xmax, length.out=100))
>
> predicted\$heightIn <- predict(model, predicted)
> sp <- ggplot(heightweight, aes(x=ageYear, y=heightIn)) +
+     geom_point(colour="grey40")
>
> sp + geom_line(data=predicted, size=1)
>
> predictvals <- function(model, xvar, yvar, xrange=NULL, samples=100, ...) {
+     if (is.null(xrange)) {
+         if (any(class(model) %in% c("lm", "glm")))
+             xrange <- range(model\$model[[xvar]])
+         else if (any(class(model) %in% "loess"))
+             xrange <- range(model\$x)
+     }
+
+ newdata <- data.frame(x = seq(xrange, xrange, length.out = samples))
+ names(newdata) <- xvar
+ newdata[[yvar]] <- predict(model, newdata = newdata, ...)
+ newdata
+ }
> modlinear <- lm(heightIn ~ ageYear, heightweight)
> modloess <- loess(heightIn ~ ageYear, heightweight)
>
> lm_predicted <- predictvals(modlinear, "ageYear", "heightIn")
> loess_predicted <- predictvals(modloess, "ageYear", "heightIn")
> sp + geom_line(data=lm_predicted, colour="red", size=.8) +
+     geom_line(data=loess_predicted, colour="blue", size=.8)
> md <- data.frame(deaths = as.numeric(mdeaths),
+                  month = as.numeric(cycle(mdeaths)))
> md <- ddply(md, "month", summarise, deaths = mean(deaths))
> p <- ggplot(md, aes(x=month, y=deaths)) + geom_line() +
+     scale_x_continuous(breaks=1:12)

> p + coord_polar()
>
> p + coord_polar() + ylim(0, max(md\$deaths))
> mdx <- md[md\$month==12, ]
>
> mdx\$month <- 0
> mdnew <- rbind(mdx, md)
>
> p %+% mdnew + coord_polar() + ylim(0, max(md\$deaths))
>
> str(economics)
Classes ‘tbl_df’, ‘tbl’ and 'data.frame': 574 obs. of  6 variables:
\$ date    : Date, format:  ...
\$ pce     : num  507 510 516 513 518 ...
\$ pop     : int  198712 198911 199113 199311 199498 199657 199808 199920 200056 200208 ...
\$ psavert : num  12.5 12.5 11.7 12.5 12.5 12.1 11.7 12.2 11.6 12.2 ...
\$ uempmed : num  4.5 4.7 4.6 4.9 4.7 4.8 5.1 4.5 4.1 4.6 ...
\$ unemploy: int  2944 2945 2958 3143 3066 3018 2878 3001 2877 2709 ...
> ggplot(economics, aes(x=date, y=psavert)) + geom_line()
>
> econ <- subset(economics, date >= as.Date("1992-05-01") &
+                    date < as.Date("1993-06-01"))
> p <- ggplot(econ, aes(x=date, y=psavert)) + geom_line()
> datebreaks <- seq(as.Date("1992-06-01"), as.Date("1993-06-01"), by="2 month")
>
> p + scale_x_date(breaks=datebreaks) +
+     theme(axis.text.x = element_text(angle=30, hjust=1))
> library(scales)
Warning message:
package ‘scales’ was built under R version 3.5.2
> library(scales)
> p + scale_x_date(breaks=datebreaks, labels=date_format("%Y %b")) +
+     theme(axis.text.x = element_text(angle=30, hjust=1))
>
> www <- data.frame(minute = as.numeric(time(WWWusage)),
+                   users = as.numeric(WWWusage))
> timeHM_formatter <- function(x) {
+     h <- floor(x/60)
+     m <- floor(x %% 60)
+     lab <- sprintf("%d:%02d", h, m) # Format the strings as HH:MM
+     return(lab)
+ }
>
> ggplot(www, aes(x=minute, y=users)) + geom_line()
> ggplot(www, aes(x=minute, y=users)) + geom_line() +
+     scale_x_continuous(name="time", breaks=seq(0, 100, by=10),
+                        labels=timeHM_formatter)
> scale_x_continuous(breaks=c(0, 20, 40, 60, 80, 100),
+                    labels=c("0:00", "0:20", "0:40", "1:00", "1:20", "1:40"))

<ScaleContinuousPosition>
Range:
Limits:    0 --    1
> timeHM_formatter(c(0, 50, 51, 59, 60, 130, 604))
 "0:00"  "0:50"  "0:51"  "0:59"  "1:00"  "2:10"
 "10:04"
>
> timeHMS_formatter <- function(x) {
+     h <- floor(x/3600)
+     m <- floor((x/60) %% 60)
+     s <- round(x %% 60) # Round to nearest second
+     lab <- sprintf("%02d:%02d:%02d", h, m, s)
+ lab <- sub("^00:", "", lab) # Remove leading 00: if present
+ lab <- sub("^0", "", lab) # Remove leading 0 if present
+ return(lab)
+ }
> islice <- subset(isabel, z == min(z))
>
> ggplot(islice, aes(x=x, y=y)) +
+     geom_segment(aes(xend = x + vx/50, yend = y + vy/50),
+                  size = 0.25)

### Latin Cubes

What is Latin Cubes? The practical applications of Latin crops and related designs are factorial experiment.factorial experiments for me... 