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)
[1] 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)[1]
> 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)
[1] 9
> f(a=1,b=2,c=3,drop=FALSE)
mvp object algebraically equal to
9
> x == mpoly_to_mvp(mpoly::as.mpoly(x))
[1] 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)[1]
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")
Loading required package: timeDate
Loading required package: timeSeries
Loading required package: fBasics
Loading required package: fAssets
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,  :
  object 'del0' not found
> rdonlp2NLP(start, objective,
+            lower=0, upper=1, linCons, funCons, control=list())
Error in rdonlp2NLP(start, objective, lower = 0, upper = 1, linCons, funCons,  :
  object 'linCons' not found
> 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`
[1] 500

$iter.max
[1] 400

$trace
[1] 0

$abs.tol
[1] 1e-20

$rel.tol
[1] 1e-10

$x.tol
[1] 1.5e-08

$step.min
[1] 2.2e-14

$scale
[1] 1

$R
[1] 1

$beta.tol
[1] 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()   
[1] NA
> amplNLPControl(
+     solver="minos", project="ampl", trace=FALSE)
$`solver`
[1] "minos"

$project
[1] "ampl"

$trace
[1] FALSE

> rdonlp2NLP(start, objective,
+            lower=0, upper=1, linCons, funCons, control=list()) 
Error in rdonlp2NLP(start, objective, lower = 0, upper = 1, linCons, funCons,  :
  object 'linCons' not found
> 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,  :
  object 'del0' not found
> 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,  :
  object 'linCons' not found
> 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,  :
  object 'objective' not found
> solnpNLPControl(
+     rho=1, outer.iter=400, inner.iter=800, delta=1e-07, tol=1e-08, trace=0)
$`rho`
[1] 1

$outer.iter
[1] 400

$inner.iter
[1] 800

$delta
[1] 1e-07

$tol
[1] 1e-08

$trace
[1] 0

> rsolnp
Error: object 'rsolnp' not found
>
> rnlminb2NLP(start, objective,
+             lower=0, upper=1, linCons, funCons, control=list()) 
Error in rnlminb2NLP(start, objective, lower = 0, upper = 1, linCons,  :
  object 'linCons' not found
> 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`
[1] 500

$iter.max
[1] 400

$trace
[1] 0

$abs.tol
[1] 1e-20

$rel.tol
[1] 1e-10

$x.tol
[1] 1.5e-08

$step.min
[1] 2.2e-14

$scale
[1] 1

$R
[1] 1

$beta.tol
[1] 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()   
[1] NA
> amplNLPControl(
+     solver="minos", project="ampl", trace=FALSE)
$`solver`
[1] "minos"

$project
[1] "ampl"

$trace
[1] FALSE

> rdonlp2NLP(start, objective,
+            lower=0, upper=1, linCons, funCons, control=list()) 
Error in rdonlp2NLP(start, objective, lower = 0, upper = 1, linCons, funCons,  :
  object 'linCons' not found
> 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,  :
  object 'del0' not found
> 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`
[1] 500

$iter.max
[1] 400

$trace
[1] 0

$abs.tol
[1] 1e-20

$rel.tol
[1] 1e-10

$x.tol
[1] 1.5e-08

$step.min
[1] 2.2e-14

$scale
[1] 1

$R
[1] 1

$beta.tol
[1] 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) :
  object 'object' not found
>
> backtestAssetsPlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestAssetsPlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
  object 'object' not found
> backtestWeightsPlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestWeightsPlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
  object 'object' not found
> backtestRebalancePlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestRebalancePlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
  object 'object' not found
> backtestPortfolioPlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestPortfolioPlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
  object 'object' not found
> backtestDrawdownPlot(object, labels=TRUE, legend=TRUE, at=NULL, format=NULL)
Error in backtestDrawdownPlot(object, labels = TRUE, legend = TRUE, at = NULL,  :
  object 'object' not found
> backtestReportPlot(object, cex=0.6, font=1, family="mono")
Error in getStrategyFun(object$backtest) : object 'object' not found
> backtestStats(100, FUN = "rollingSigma", ...)
Error: '...' used in an incorrect context
> library("PortfolioAnalytics", lib.loc="~/R/win-library/3.5")
Loading required package: zoo

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

Loading required package: xts
Loading required package: foreach
Loading required package: PerformanceAnalytics

Package PerformanceAnalytics (1.5.2) loaded.
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)
> head(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)
[1] "Zurich"
> finCenter(NYC)
[1] "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")
> pspec <- add.constraint(pspec, type="box", min=0.05, max=0.45)
> pspec <- add.constraint(pspec,
+                         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) :
  object 'buckets' not found
> 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
[1] NA
> library("FRAPO", lib.loc="~/R/win-library/3.5")
Loading required package: cccp
Loading required package: Rglpk
Loading required package: slam
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)
> head(cs)
[1] 200 200 200 200 200 200
> yret <- diff(log(y))
> bilson <- trdbilson(yret, exponent = 2)
> head(bilson)
[1]  0.01948730 -0.01935785  0.01177193 -0.04517835
[5]  0.10349955 -0.01545405
> es <- trdes(yret, lambda = 0.95)
> head(es)
[1]  0.018485339 -0.017438597  0.010304526
[4] -0.042136563  0.093941299 -0.009969699
> cs <- capser(y, min = 100, max = 200)
> yret <- diff(log(y))
> binary <- trdbinary(yret)
> head(binary)
[1]  0.02477189 -0.02460780  0.01497858 -0.05712579
[5]  0.12829276 -0.01965560
> hp <- trdhp(y, lambda = 1600)
> head(hp)
[1] 387.1268 389.8778 392.6293 395.3850 398.1460
[6] 400.9114
> sma <- trdsma(y, n.periods = 24)
> head(sma, 30)
 [1]       NA       NA       NA       NA       NA
 [6]       NA       NA       NA       NA       NA
[11]       NA       NA       NA       NA       NA
[16]       NA       NA       NA       NA       NA
[21]       NA       NA       NA 418.2483 420.7617
[26] 423.6004 426.5617 429.7025 433.3096 435.5387
> wma <- trdwma(y, weights = c(0.4, 0.3, 0.2, 0.1))
> head(wma, 30)
 [1]      NA      NA      NA 391.205 384.938
 [6] 395.898 402.343 408.264 408.556 410.505
[11] 412.633 411.220 416.691 416.362 417.228
[16] 418.201 423.188 430.449 435.296 439.885
[21] 445.320 444.131 446.805 448.474 448.468
[26] 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))
[1] 50 25  0  0  0  0 25
> zapsmall(convolve(x, y[3:5], type = "f"))
[1]  0 25 50 25  0
> x <- rnorm(50)
> y <- rnorm(50)
> all.equal(convolve(x, y, conj = FALSE), rev(convolve(rev(y),x)))
[1] 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
[1] 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] 1.033928 1.227142 1.381714 1.505371 1.604297
> sqrt(a.pred$h)
[1] 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[1]> plot(x, y)
Browse[1]>
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")
Error in xy.coords(x, y) : object 'yS' not found
> 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] 1

$alpha
[1] 0.3

$beta
[1] 0.5

$esp0
[1] -1.5

$h0
[1] 4.96

attr(,"class")
[1] "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
[1] 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")
[1] 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
[1] 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
[1] 3.287939e-09
> abs((res2a - res2)) # 5.331195e-11, intercept/slope better numerically
[1] 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
[1] 6.938894e-18
> abs((res2b - res2)) # 1.040834e-16
[1] 1.040834e-16
> abs((res1b - res2)/res2) # 2.6119e-16 # both within machine precision
[1] 2.6119e-16
> abs((res2b - res2)/res2) # 3.91785e-15
[1] 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[1], sigma = x[2], weights = 1)
+ max(abs((vPA - vAz1))) # 5.551115e-17
Error: unexpected symbol in:
"                              mu = x[1], sigma = x[2], weights = 1)
max"
> max(abs((vPA - vAz2))) #   ""
Error: object 'vPA' not found
>
> max(abs((vPA - vAz1a))) # 3.287941e-09
Error: object 'vPA' not found
> max(abs((vPA - vAz2a))) #  1.465251e-10, intercept/slope better
Error: object 'vPA' not found
>
> max(abs((vPA - vAz1b))) # 4.374869e-13
Error: object 'vPA' not found
> max(abs((vPA - vAz2b))) # 3.330669e-16
Error: object 'vPA' not found
> }
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",  :
  object 'hw' not found
>
> 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) +
+     scale_fill_gradient(low="black", high="white")
>
> ggplot(heightweight, aes(x=weightLb, y=heightIn, fill=ageYear)) +
+     geom_point(shape=21, size=2.5) +
+     scale_fill_gradient(low="black", high="white", breaks=12:17,
+                         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) +
+     scale_fill_gradient(low="lightblue", high="red", limits=c(0, 6000))
>
> library(hexbin)
Warning message:
package ‘hexbin’ was built under R version 3.5.2
> library(hexbin)
> sp + stat_binhex() +
+     scale_fill_gradient(low="lightblue", high="red",
+                         limits=c(0, 8000))
>
> sp + stat_binhex() +
+     scale_fill_gradient(low="lightblue", high="red",
+                         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[1], xrange[2], 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))
[1] "0:00"  "0:50"  "0:51"  "0:59"  "1:00"  "2:10"
[7] "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...