Monday, 29 April 2019

FX derivatives can be grouped in three main categories

FX derivatives can be grouped in three main categories:

Futures, Swaps, and Options.  we will only deal with option-type derivatives.
The price of the underlying will always be S, but if it is not necessarily a
currency, we will use numeric or alphabetic indices such as S1
or SA
• The strike price will always be X
• The expected value operator will be denoted by E

Currency options

European currency options grant the holder the right to buy (call option) or sell
(put option) currency at a predetermined exchange rate (strike price or exercise price, X), on a specified date (maturity, T). These financial assets are also called foreign
exchange options (or FX options),
We can use the Black-Scholes option or the GBSOption function, both of which are practically the same, and the latter is a shorthand alias
for the prior function.
BlackScholesOption(TypeFlag, S, X, Time, r, b, sigma,...)
Here, TypeFlag is either the character c, which stands for call or p (put). S is the
current price and sigma is the volatility of the underlying. X is the strike price and
Time is the time to maturity
The other two parameters are a bit tricky because r and b are the risk-free rates
we must set b = r to get the BS stock option model,
and set b = r-q to get the currency option model or the stock option model with
continuous dividend yield
GBSOption(TypeFlag = "c", S = 0.745, X = 0.7, Time = 5, r = 0.03, b = 0.01, sigma = 0.2)
BlackScholesOption ("c", 0.7450, 0.7, 5, 0.03, 0.01, 0.2)

Title:
 Black Scholes Option Valuation

Call:
 GBSOption(TypeFlag = "c", S = 0.745, X = 0.7, Time = 5, r = 0.03,
     b = 0.01, sigma = 0.2)

Parameters:
          Value:
 TypeFlag c   
 S        0.745
 X        0.7 
 Time     5   
 r        0.03
 b        0.01
 sigma    0.2 

Option Price:
 0.152222

Description:
 Mon Apr 29 16:34:31 2019
We can also check out the price of the put option:
 BlackScholesOption("p", 0.7450, 0.7, 5, 0.03, 0.01, 0.2)

Title:
 Black Scholes Option Valuation 

Call:
 GBSOption(TypeFlag = "p", S = 0.745, X = 0.7, Time = 5, r = 0.03, 
     b = 0.01, sigma = 0.2)

Parameters:
          Value:
 TypeFlag p     
 S        0.745 
 X        0.7   
 Time     5     
 r        0.03  
 b        0.01  
 sigma    0.2   

Option Price:
 0.08061367 

Description:

 Mon Apr 29 16:35:32 2019 
we can also check the consistency with the put-call parity
> c - p = S*exp(-r*T)–X*exp(-q*T)
> c - p = 0.152222 - 0.08061367=0.07160833
On the right-hand side, we have: 0.745*exp(-0.02*5)-0.7*exp(-0.03*5) = 0.07160829
Correlation between the Wiener processes changes the picture dramatically. In the
case of positive correlation, the two Wiener processes look like they are moving in
the same direction; in the case of negative correlation, they look like they are moving

in the opposite direction.
D2_Wiener <- function() {
+     dev.new(width = 10, height = 4)
+     par(mfrow = c(1, 3), oma = c(0, 0, 2, 0))
+     for(i in 1:3) {
+         W1 <- cumsum(rnorm(100000))
+         W2 <- cumsum(rnorm(100000))
+         plot(W1,W2, type= "l", ylab = "", xlab = "")
+     }
+     mtext("2-dimensional Wiener-processes with no correlation",
+           outer = TRUE, cex = 1.5, line = -1)
+ }

https://www.mathclasstutor.online

https://www.mathclasstutor.online

> D2_Wiener()
> Correlated_Wiener <- function(cor) {
+     dev.new(width = 10, height = 4)
+     par(mfrow = c(1, 3), oma = c(0, 0, 2, 0))
+     for(i in 1:3) {
+         W1 <- cumsum(rnorm(100000))
+         W2 <- cumsum(rnorm(100000))
+         W3 <- cor * W1 + sqrt(1 - cor^2) * W2
+         plot(W1, W3, type= "l", ylab = "", xlab = "")
+     }
+     mtext(paste("2-dimensional Wiener-processes (",cor," correlation)",
+                 sep = ""), outer = TRUE, cex = 1.5, line = -1)
+ }
> Correlated_Wiener(0.6)

> Correlated_Wiener(-0.7)
Margrabe formula.
Margrabe <- function(S1, S2, sigma1, sigma2, Time, rho, delta1 = 0,
+                      delta2 = 0) {
+     sigma <- sqrt(sigma1^2 + sigma2^2 - 2 * sigma1 * sigma2 * rho)
+     d1 <- ( log(S1/S2) + ( delta2-delta1 + sigma^2/2 ) * Time ) /
+         (sigma*sqrt(Time))
+     d2 <- ( log(S1/S2) + ( delta2-delta1 - sigma^2/2 ) * Time ) /
+         (sigma*sqrt(Time))
+     M <- S1*exp(-delta1*Time)*pnorm(d1) - S2*exp(-delta2*Time)*pnorm(d2)
+     return(M)

+ }
Margrabe(100, 120, .2, .3, 2, .15)
[1] 12.05247
> Margrabe(100, 120, .2, 0, 2, 0, 0, 0.03)

[1] 6.566047
 we calculate the Margrabe price of the option for different values of correlation.
https://www.mathclasstutor.online

The result is not surprising. When the correlation is high, we have the right to switch
between identical stocks, which, clearly, is worth nothing. When the correlation
is high on the negative side, we have better chances to make a good deal with the
option if things go wrong (which means that if our asset decreases, the higher the negative correlation, the higher the chance that the price of the other asset increases and saves us from loss). In other words, in this case, the option is for insurance rather
than speculation; we do not have to bear the risk from the price change of the other
asset. This is why the option is more valuable when the correlation is negative

Thursday, 25 April 2019

PerformanceAnalytics

PerformanceAnalytics

library(PerformanceAnalytics)
> library ( quadprog )
> library ( rrcov )
Loading required package: robustbase
Scalable Robust Estimators with High Breakdown Point (version 1.4-7)

Warning messages:
1: package ‘rrcov’ was built under R version 3.5.2
2: package ‘robustbase’ was built under R version 3.5.1
> library ( rrcov )
> Dmat       <- matrix(0,3,3)
> diag(Dmat) <- 1
> dvec       <- c(0,5,0)
> Amat       <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
> bvec       <- c(-8,2,0)
> solve.QP(Dmat,dvec,Amat,bvec=bvec)
$`solution`
[1] 0.4761905 1.0476190 2.0952381

$value
[1] -2.380952

$unconstrained.solution
[1] 0 5 0

$iterations
[1] 3 0

$Lagrangian
[1] 0.0000000 0.2380952 2.0952381

$iact
[1] 3 2
https://www.mathclasstutor.online


> library ( zoo )
> data ( StockIndex )
> pzoo <− zoo ( StockIndex , order.by = rownames( StockIndex ) )
> rzoo <− (pzoo/lag ( pzoo , k = −1) − 1 ) ∗ 100
> boxplot ( coredata ( rzoo ) )
> rstats <− rbind ( apply ( rzoo , 2 , summary ) ,
+                   skewness ( rzoo ) ,
+                   kurtosis ( rzoo )
+ )

Sunday, 21 April 2019

The functions are exported via the NAMESPACE mechanism & Gumbel Hougaard Copula

The functions are exported via the NAMESPACE mechanism

You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Workspace loaded from ~/.RData]

Loading required package: fOptions
Loading required package: timeDate
Loading required package: timeSeries
Loading required package: fBasics
> install.packages("BLCOP")
Installing package into ‘C:/Users/ADMIN/Documents/R/win-library/3.5’
(as ‘lib’ is unspecified)
also installing the dependency ‘RUnit’

trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.5/RUnit_0.4.32.zip'
Content type 'application/zip' length 296247 bytes (289 KB)
downloaded 289 KB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.5/BLCOP_0.3.1.zip'
Content type 'application/zip' length 811518 bytes (792 KB)
downloaded 792 KB

package ‘RUnit’ successfully unpacked and MD5 sums checked
package ‘BLCOP’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
C:\Users\ADMIN\AppData\Local\Temp\RtmpM3bcRt\downloaded_packages
Warning messages:
1: package ‘fOptions’ 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.2
> library(BLCOP)
Loading required package: MASS
Loading required package: quadprog
Warning message:
package ‘BLCOP’ was built under R version 3.5.3
> pick <- newPMatrix(letters[1:8], 3)
> pick[1,7] <- 1
> pick[2,1] <- -1
> pick[2,2] <- 1
> pick[3, 3:6] <- c(0.9, -0.9, .1, -.1)
> confidences <- 1 / c(0.00709, 0.000141, 0.000866)
> myViews <- BLViews(pick, q = c(0.0525, 0.0025, 0.02), confidences, letters[1:8]
+ )
> myViews
1 : 1*g=0.0525  + eps. Confidence: 141.0437
2 : -1*a+1*b=0.0025  + eps. Confidence: 7092.199
3 : 0.9*c+-0.9*d+0.1*e+-0.1*f=0.02  + eps. Confidence: 1154.734
> dispersion <- c(.376,.253,.360,.333,.360,.600,.397,.396,.578,.775) / 1000
> sigma <- BLCOP:::.symmetricMatrix(dispersion, dim = 4)
> caps <- rep(1/4, 4)
> mu <- 2.5 * sigma
> dim(mu) <- NULL
> marketDistribution <- mvdistribution("mt", mean = mu, S = sigma, df = 5 )
Error in mvdistribution("mt", mean = mu, S = sigma, df = 5) :
  Sampling function for this distribution does not exist!
> pick <- newPMatrix(c("SP", "FTSE", "CAC", "DAX"), 1)
> pick[1,4] <- 1
> vdist <- list(distribution("unif", min = -0.02, max = 0))
> views <- COPViews(pick, vdist, 0.2, c("SP", "FTSE", "CAC", "DAX"))
> CAPMList(monthlyReturns, marketIndex = sp500Returns, riskFree = US13wTB, regFunc = "rlm")
          alphas     betas
IBM  0.015540244 1.1901296
MS   0.053239036 1.7398334
DELL 0.023938008 1.3389280
C    0.016934353 1.1822813
JPM  0.035134257 1.5712330
BAC  0.007196188 0.9702085
> x <- c(0.001005,0.001328,-0.000579,-0.000675,0.000121,0.000128,-0.000445,-0.000437 ,
+        0.001328,0.007277,-0.001307,-0.000610,-0.002237,-0.000989,0.001442,-0.001535 ,
+        -0.000579,-0.001307,0.059852,0.027588,0.063497,0.023036,0.032967,0.048039 ,
+        -0.000675,-0.000610,0.027588,0.029609,0.026572,0.021465,0.020697,0.029854 ,
+        0.000121,-0.002237,0.063497,0.026572,0.102488,0.042744,0.039943,0.065994 ,
+        0.000128,-0.000989,0.023036,0.021465,0.042744,0.032056,0.019881,0.032235 ,
+        -0.000445,0.001442,0.032967,0.020697,0.039943,0.019881,0.028355,0.035064 ,
+        -0.000437,-0.001535,0.048039,0.029854,0.065994,0.032235,0.035064,0.079958 )
> varCov <- matrix(x, ncol = 8, nrow = 8)
> mu <- c(0.08, 0.67,6.41, 4.08, 7.43, 3.70, 4.80, 6.60) / 100
> pick <- matrix(0, ncol = 8, nrow = 3, dimnames = list(NULL, letters[1:8]))
> pick[1,7] <- 1
> pick[2,1] <- -1; pick[2,2] <- 1
> pick[3, 3:6] <- c(0.9, -0.9, .1, -.1)
> confidences <- 1 / c(0.000709, 0.000141, 0.000866)
> myViews <- BLViews(pick, c(0.0525, 0.0025, 0.02), confidences, letters[1:8])
> myPosterior <- posteriorEst(myViews, tau = 0.025, mu, varCov )
> myPosterior
Prior means:
     a      b      c      d      e
0.0008 0.0067 0.0641 0.0408 0.0743
     f      g      h
0.0370 0.0480 0.0660
Posterior means:
          a           b           c
0.000636636 0.004983111 0.064945904
          d           e           f
0.043209911 0.075889579 0.039353875
          g           h
0.049335088 0.068439660
Posterior covariance:
              [,1]          [,2]
[1,]  0.0010297637  0.0013570700
[2,]  0.0013570700  0.0073802381
[3,] -0.0005877981 -0.0013326418
[4,] -0.0006868193 -0.0006311548
[5,]  0.0001321052 -0.0022651315
[6,]  0.0001366368 -0.0010033018
[7,] -0.0004518421  0.0014486393
[8,] -0.0004403174 -0.0015640550
              [,3]          [,4]
[1,] -0.0005877981 -0.0006868193
[2,] -0.0013326418 -0.0006311548
[3,]  0.0606199378  0.0280331207
[4,]  0.0280331207  0.0301452738
[5,]  0.0642059208  0.0269354857
[6,]  0.0233376823  0.0218124452
[7,]  0.0333195421  0.0209724486
[8,]  0.0486120217  0.0303038126
              [,5]          [,6]
[1,]  0.0001321052  0.0001366368
[2,] -0.0022651315 -0.0010033018
[3,]  0.0642059208  0.0233376823
[4,]  0.0269354857  0.0218124452
[5,]  0.1039853518  0.0434735925
[6,]  0.0434735925  0.0326741247
[7,]  0.0403758705  0.0201390307
[8,]  0.0668800713  0.0327345035
              [,7]          [,8]
[1,] -0.0004518421 -0.0004403174
[2,]  0.0014486393 -0.0015640550
[3,]  0.0333195421  0.0486120217
[4,]  0.0209724486  0.0303038126
[5,]  0.0403758705  0.0668800713
[6,]  0.0201390307  0.0327345035
[7,]  0.0286908498  0.0354784746
[8,]  0.0354784746  0.0813542571
> confidences(myViews) <-  1 / c(0.0001, 0.0001, 0.0005)
> myPosterior2 <- posteriorEst(myViews, tau = 0.025, mu, varCov )
> myPosterior2
Prior means:
     a      b      c      d      e
0.0008 0.0067 0.0641 0.0408 0.0743
     f      g      h
0.0370 0.0480 0.0660
Posterior means:
           a            b            c
0.0005567763 0.0046938861 0.0668960593
           d            e            f
0.0455483701 0.0784192403 0.0415007851
           g            h
0.0516138866 0.0713497046
Posterior covariance:
              [,1]          [,2]
[1,]  0.0010295987  0.0013563765
[2,]  0.0013563765  0.0073668897
[3,] -0.0005828989 -0.0013321345
[4,] -0.0006820877 -0.0006324597
[5,]  0.0001383777 -0.0022612469
[6,]  0.0001410988 -0.0010018496
[7,] -0.0004468617  0.0014431554
[8,] -0.0004338236 -0.0015631163
              [,3]          [,4]
[1,] -0.0005828989 -0.0006820877
[2,] -0.0013321345 -0.0006324597
[3,]  0.0602995605  0.0278395293
[4,]  0.0278395293  0.0299654279
[5,]  0.0638145738  0.0266953375
[6,]  0.0231479118  0.0216471782
[7,]  0.0330542399  0.0207699549
[8,]  0.0482749436  0.0300556089
              [,5]          [,6]
[1,]  0.0001383777  0.0001410988
[2,] -0.0022612469 -0.0010018496
[3,]  0.0638145738  0.0231479118
[4,]  0.0266953375  0.0216471782
[5,]  0.1035062513  0.0432382103
[6,]  0.0432382103  0.0325205456
[7,]  0.0400505520  0.0199488216
[8,]  0.0664653327  0.0324987249
              [,7]          [,8]
[1,] -0.0004468617 -0.0004338236
[2,]  0.0014431554 -0.0015631163
[3,]  0.0330542399  0.0482749436
[4,]  0.0207699549  0.0300556089
[5,]  0.0400505520  0.0664653327
[6,]  0.0199488216  0.0324987249
[7,]  0.0284411037  0.0351700306
[8,]  0.0351700306  0.0809678994
> dispersion <- c(.376,.253,.360,.333,.360,.600,.397,.396,.578,.775) / 1000
> sigma <- BLCOP:::.symmetricMatrix(dispersion, dim = 4)
> caps <- rep(1/4, 4)
> mu <- 2.5 * sigma
> dim(mu) <- NULL
> marketDistribution <- mvdistribution("mt", mean = mu, S = sigma, df = 5 )
Error in mvdistribution("mt", mean = mu, S = sigma, df = 5) :
  Sampling function for this distribution does not exist!
> pick <- matrix(0, ncol = 4, nrow = 1, dimnames = list(NULL, c("SP", "FTSE", "CAC", "DAX"))
+ )
> pick[1,4] <- 1
> vdist <- list(distribution("unif", min = -0.02, max = 0))
> views <- COPViews(pick, vdist, 0.2, c("SP", "FTSE", "CAC", "DAX"))
> posterior <- COPPosterior(marketDistribution, views)
Error in classNames %in% class(object) :
  object 'marketDistribution' not found
> stocks <- colnames(monthlyReturns)
> pick <- matrix(0, ncol = 6, nrow = 2, dimnames = list(NULL, stocks))
> pick[1,"IBM"] <- 1
> pick[1, "DELL"] <- 0.04
> pick[2, "C"] <- 1
> pick[2, "JPM"] <- 0.6
> confidences <- 1 / c(0.7, 0.1)
> views <- BLViews( P = pick, q = c(0.1,0.1) , confidences = confidences,stocks)
> deleteViews(views, 1)
1 : 1*C+0.6*JPM=0.1  + eps. Confidence: 10
> dispersion <- c(.376,.253,.360,.333,.360,.600,.397,.396,.578,.775) / 1000
> sigma <- BLCOP:::.symmetricMatrix(dispersion, dim = 4)
> caps <- rep(1/4, 4)
> mu <- 2.5 * sigma
> dim(mu) <- NULL
> marketDistribution <- mvdistribution("mt", mean = mu, S = sigma, df = 5 )
Error in mvdistribution("mt", mean = mu, S = sigma, df = 5) :
  Sampling function for this distribution does not exist!
> pick <- matrix(0, ncol = 4, nrow = 1, dimnames = list(NULL, c("SP", "FTSE", "CAC", "DAX")))
> pick[1,4] <- 1
> vdist <- list(distribution("unif", min = -0.02, max = 0))
> views <- COPViews(pick, vdist, 0.2, c("SP", "FTSE", "CAC", "DAX"))
> posterior <- COPPosterior(marketDistribution, views)
Error in classNames %in% class(object) :
  object 'marketDistribution' not found
> densityPlots(posterior, 4)
Error in densityPlots(posterior, 4) : object 'posterior' not found
> ts.plot(sp500Returns)
> ts.plot(US13wTB)
> x <- distribution("pois", lambda = 5)
> hist(sampleFrom(x, 1000), col = "blue", prob = TRUE)
> myUnif <- distribution("unif", min = -0.1, max = 0.1)
> hist(sampleFrom(myUnif, 1000))
> mvNormal <- mvdistribution("mnorm", mean = c(1, 5), varcov = diag(c(2, 0.1)))
Error in mvdistribution("mnorm", mean = c(1, 5), varcov = diag(c(2, 0.1))) :
  Sampling function for this distribution does not exist!
> x <- sampleFrom(mvNormal, 1000)
Error in classNames %in% class(object) : object 'mvNormal' not found
> mvNormal <- mvdistribution("mnorm", mean = c(1, 5), varcov = diag(c(2, 0.1)))
Error in mvdistribution("mnorm", mean = c(1, 5), varcov = diag(c(2, 0.1))) :
  Sampling function for this distribution does not exist!
> showClass("distribution")
Class "distribution" [package "BLCOP"]

Slots:
                         
Name:       RName parameters

Class:  character    numeric

> showClass("mvdistribution")
Class "mvdistribution" [package "BLCOP"]

Slots:
                         
Name:       RName parameters
Class:  character       list
> showClass("distribution")
Class "distribution" [package "BLCOP"]

Slots:
                         
Name:       RName parameters
Class:  character    numeric
> CAPMList(monthlyReturns, marketIndex = sp500Returns, riskFree = US13wTB)
          alphas     betas
IBM  0.020883598 1.2806423
MS   0.059548398 1.8186393
DELL 0.017010062 1.2843338
C    0.014492325 1.1953415
JPM  0.027365230 1.5623610
BAC  0.002829908 0.9628988
> pickMatrix <- matrix(c(rep(1/2, 2), -1,  rep(0, 3)), nrow = 1, ncol = 6 )
> views <- BLViews(P = pickMatrix, q = 0.08,confidences =  100,
+                  assetNames = colnames(monthlyReturns))
> marketPosterior <- BLPosterior(monthlyReturns, views, marketIndex = sp500Returns,
+                                riskFree = US13wTB)
> posteriorFeasibility(marketPosterior)
$`mahalDist`
[1] 0.06484571

$mahalDistProb
[1] 0.9999945

$sensitivities
[1] "Not implemented yet"

> install.packages("gumbel")
Installing package into ‘C:/Users/ADMIN/Documents/R/win-library/3.5’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.5/gumbel_1.10-2.zip'
Content type 'application/zip' length 471003 bytes (459 KB)
downloaded 459 KB

package ‘gumbel’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
C:\Users\ADMIN\AppData\Local\Temp\RtmpM3bcRt\downloaded_packages
> library(gumbel)
Warning message:
package ‘gumbel’ was built under R version 3.5.2
> u<-seq(0,1, .1)
> v<-u
> dgumbel(u,v,1)
 [1] NaN   1   1   1   1   1   1   1   1
[10]   1 NaN
> pgumbel(u,v,1)
 [1] 0.00 0.01 0.04 0.09 0.16 0.25 0.36
 [8] 0.49 0.64 0.81 1.00
> dgumbel(cbind(u,v), alpha=1)
 [1] NaN   1   1   1   1   1   1   1   1
[10]   1 NaN
> pgumbel(cbind(u,v), alpha=1)
 [1] 0.00 0.01 0.04 0.09 0.16 0.25 0.36
 [8] 0.49 0.64 0.81 1.00
> x <- cbind(u,u,u)
> y <- array(u, c(1,11,3))
> pgumbel(x, alpha=2, dim=3)
 [1] 0.00000000 0.01853315 0.06156706
 [4] 0.12426461 0.20452561 0.30102374
 [7] 0.41280666 0.53914047 0.67943346
[10] 0.83319317 1.00000000
> pgumbel(y, alpha=2, dim=3)
     [,1]       [,2]       [,3]
[1,]    0 0.01853315 0.06156706
          [,4]      [,5]      [,6]
[1,] 0.1242646 0.2045256 0.3010237
          [,7]      [,8]      [,9]
[1,] 0.4128067 0.5391405 0.6794335
         [,10] [,11]
[1,] 0.8331932     1
> dgumbel(x, alpha=2, dim=3)
 [1]       NaN  6.922376  3.646759
 [4]  2.770993  2.510089  2.586175
 [7]  3.024395  4.149303  7.366620
[10] 23.650259       NaN
> dgumbel(y, alpha=2, dim=3)
     [,1]     [,2]     [,3]     [,4]
[1,]  NaN 6.922376 3.646759 2.770993
         [,5]     [,6]     [,7]     [,8]
[1,] 2.510089 2.586175 3.024395 4.149303
        [,9]    [,10] [,11]
[1,] 7.36662 23.65026   NaN
> x <- cbind(x,u)
> pgumbel(x, alpha=3, dim=4)
 [1] 0.00000000 0.02585824 0.07770595
 [4] 0.14790462 0.23351222 0.33277038
 [7] 0.44446448 0.56768637 0.70172176
[10] 0.84598860 1.00000000
> y <- array(u, c(2,1,11,4))
> pgumbel(y, alpha=3, dim=4)
, , 1

           [,1]
[1,] 0.00000000
[2,] 0.02585824

, , 2

           [,1]
[1,] 0.07770595
[2,] 0.14790462

, , 3

          [,1]
[1,] 0.2335122
[2,] 0.3327704

, , 4

          [,1]
[1,] 0.4444645
[2,] 0.5676864

, , 5

          [,1]
[1,] 0.7017218
[2,] 0.8459886

, , 6

     [,1]
[1,]    1
[2,]    0

, , 7

           [,1]
[1,] 0.02585824
[2,] 0.07770595

, , 8

          [,1]
[1,] 0.1479046
[2,] 0.2335122

, , 9

          [,1]
[1,] 0.3327704
[2,] 0.4444645

, , 10

          [,1]
[1,] 0.5676864
[2,] 0.7017218

, , 11

          [,1]
[1,] 0.8459886
[2,] 1.0000000

> rand <- t(rgumbel(200,1))
> plot(rand[1,], rand[2,], col="green", main="Gumbel copula")
> rand <- t(rgumbel(200,2))
> plot(rand[1,], rand[2,], col="red", main="Gumbel copula")
> nbsimu <- 10000
> system.time(rgumbel(nbsimu, 2, dim=2, method=1))[3]
elapsed
   0.01
> system.time(rgumbel(nbsimu, 2, dim=2, method=2))[3]
elapsed
   0.61
> anim <-function(n, max=50)
+ {
+     for(i in seq(1,max,length.out=n))
+     {
+         u <- t(rgumbel(10000, i, method=2))
+         plot(u[1,], u[2,], col="green", main="Gumbel copula",
+              xlim=c(0,1), ylim=c(0,1), pch=".")
+         cat()
+     }
+ }
> anim(20, 20)
> x <- seq(.05, .95, length = 30)
>
> y <- x
> z <- outer(x, y, dgumbel, alpha=2)
> persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
+       ltheta = 100, shade = 0.25, ticktype = "detailed",
+       xlab = "x", ylab = "y", zlab = "density")
https://www.mathclasstutor.online


> zlim <- c(0, max(z))
> ncol <- 100
> nrz <- nrow(z)
> ncz <- ncol(z)
> couleurs <- tail(jet.colors(1.2*ncol),ncol)
Error in jet.colors(1.2 * ncol) : could not find function "jet.colors"
> couleurs <- tail(colors(1.2*ncol),ncol)
> fcol <- couleurs[trunc(z/zlim[2]*(ncol-1))+1]
> dim(fcol) <- c(nrz,ncz)
> fcol <- fcol[-nrz,-ncz]
> persp(x, y, z, col=fcol, zlim=zlim, theta=30, phi=30, ticktype = "detailed",
+       box = TRUE, xlab = "x", ylab = "y", zlab = "density")
> z <- outer(x, y, pgumbel, alpha=2)
> persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
+       ltheta = 100, shade = 0.25, ticktype = "detailed",
+       xlab = "u", ylab = "v", zlab = "cdf")
> simu <- qexp(rgumbel(200, 2))
> gumbel.MBE(simu[,1], simu[,2])
[1] 1.008945 1.011470 1.951363
> gumbel.IFM(simu[,1], simu[,2])
[1] 1.014007 1.016544 2.022765
> gumbel.EML(simu[,1], simu[,2])
[1] 1.025841 1.028801 2.011885
> gumbel.CML(simu[,1], simu[,2]
+ )
[1] 1.893303
> simu <- qgamma(rgumbel(200, 3), 2, 1)
> gumbel.MBE(simu[,1], simu[,2], "gamma")
[1] 0.9459023 1.8808981 1.0406970
[4] 2.0350758 3.0133253
> gumbel.IFM(simu[,1], simu[,2], "gamma")
[1] 0.9667415 1.9224548 1.0005313
[4] 1.9565467 2.9528046
> gumbel.EML(simu[,1], simu[,2], "gamma")
[1] 0.9618522 1.9176133 1.0059094
[4] 1.9628113 2.9535937
> gumbel.CML(simu[,1], simu[,2])
[1] 2.917404

Saturday, 20 April 2019

“Finance” and “TimeSeries” fGARCH

Finance and TimeSeries

 library(bayesGARCH)
Warning message:
package ‘bayesGARCH’ was built under R version 3.5.3
> data(dem2gbp)
> y <- dem2gbp[1:750]
> MCMC <- bayesGARCH(y, control = list(n.chain = 2, l.chain = 200))
chain:  1  iteration:  10  parameters:  0.0407 0.1977 0.6747 97.9185
chain:  1  iteration:  20  parameters:  0.041 0.1566 0.6818 66.2305
chain:  1  iteration:  30  parameters:  0.0439 0.208 0.65 45.4827
chain:  1  iteration:  40  parameters:  0.0461 0.2426 0.6215 58.8274
chain:  1  iteration:  50  parameters:  0.0361 0.2158 0.6922 52.3408
chain:  1  iteration:  60  parameters:  0.0255 0.1801 0.7399 46.8184
chain:  1  iteration:  70  parameters:  0.0344 0.2125 0.7051 34.2654
chain:  1  iteration:  80  parameters:  0.037 0.1893 0.6955 35.4467
chain:  1  iteration:  90  parameters:  0.0299 0.2065 0.6975 33.0581
chain:  1  iteration:  100  parameters:  0.0332 0.2006 0.6877 28.2932
chain:  1  iteration:  110  parameters:  0.0385 0.2189 0.6584 23.1169
chain:  1  iteration:  120  parameters:  0.0377 0.1947 0.6551 18.713
chain:  1  iteration:  130  parameters:  0.0435 0.2348 0.6385 21.2452
chain:  1  iteration:  140  parameters:  0.0556 0.3015 0.5682 16.9942
chain:  1  iteration:  150  parameters:  0.0296 0.2507 0.6771 16.5932
chain:  1  iteration:  160  parameters:  0.0547 0.2298 0.5981 14.338
chain:  1  iteration:  170  parameters:  0.0502 0.2043 0.6283 10.2683
chain:  1  iteration:  180  parameters:  0.0455 0.296 0.6117 10.5093
chain:  1  iteration:  190  parameters:  0.0561 0.267 0.593 8.3902
chain:  1  iteration:  200  parameters:  0.0272 0.2558 0.6836 8.4362
chain:  2  iteration:  10  parameters:  0.0554 0.2672 0.5665 76.2903
chain:  2  iteration:  20  parameters:  0.0632 0.2456 0.5452 75.3237
chain:  2  iteration:  30  parameters:  0.0437 0.1915 0.6723 43.8814
chain:  2  iteration:  40  parameters:  0.0332 0.2959 0.6609 36.4541
chain:  2  iteration:  50  parameters:  0.0625 0.4065 0.4705 23.7808
chain:  2  iteration:  60  parameters:  0.046 0.3761 0.5178 17.951
chain:  2  iteration:  70  parameters:  0.0787 0.349 0.4474 15.7251
chain:  2  iteration:  80  parameters:  0.0671 0.2682 0.5123 20.978
chain:  2  iteration:  90  parameters:  0.0414 0.2661 0.6255 22.1921
chain:  2  iteration:  100  parameters:  0.0494 0.1933 0.6442 20.2299
chain:  2  iteration:  110  parameters:  0.0379 0.2421 0.651 23.7265
chain:  2  iteration:  120  parameters:  0.0337 0.2249 0.6749 16.5062
chain:  2  iteration:  130  parameters:  0.0288 0.2614 0.6707 15.3967
chain:  2  iteration:  140  parameters:  0.041 0.2597 0.5917 9.8975
chain:  2  iteration:  150  parameters:  0.0432 0.2261 0.6611 13.9151
chain:  2  iteration:  160  parameters:  0.031 0.2515 0.6666 12.1141
chain:  2  iteration:  170  parameters:  0.045 0.1672 0.6649 10.9084
chain:  2  iteration:  180  parameters:  0.0364 0.2269 0.7005 8.7603
chain:  2  iteration:  190  parameters:  0.044 0.2549 0.5976 10.3563
chain:  2  iteration:  200  parameters:  0.038 0.3334 0.6159 8.0582
> plot(MCMC)
> smpl <- formSmpl(MCMC, l.bi = 50)

n.chain:  2
l.chain:  200
l.bi:  50
batch.size:  1
smpl size:  300
> summary(smpl)

Iterations = 1:300
Thinning interval = 1
Number of chains = 1
Sample size per chain = 300

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean       SD Naive SE
alpha0  0.04308  0.01193 0.000689
alpha1  0.23996  0.04806 0.002774
beta    0.63368  0.06606 0.003814
nu     19.48437 10.38978 0.599854
       Time-series SE
alpha0        0.00321
alpha1        0.01294
beta          0.02203
nu            4.53866

2. Quantiles for each variable:

          2.5%      25%      50%
alpha0 0.02615  0.03506  0.04117
alpha1 0.16644  0.20486  0.23441
beta   0.46807  0.59946  0.64301
nu     8.44042 11.02220 17.59416
            75%    97.5%
alpha0  0.04805  0.07271
alpha1  0.26695  0.34679
beta    0.68334  0.72659
nu     21.86966 50.93082

> smpl <- as.matrix(smpl)
> pairs(smpl)
> MCMC <- bayesGARCH(y, lambda = 100, delta = 500,
+                    control = list(n.chain = 2, l.chain = 200))
chain:  1  iteration:  10  parameters:  0.0329 0.1778 0.6998 500.0104
chain:  1  iteration:  20  parameters:  0.0404 0.2281 0.6614 500.0087
chain:  1  iteration:  30  parameters:  0.036 0.1599 0.6875 500.0017
chain:  1  iteration:  40  parameters:  0.0545 0.1614 0.6322 500.0184
chain:  1  iteration:  50  parameters:  0.0547 0.2635 0.6081 500.0196
chain:  1  iteration:  60  parameters:  0.0466 0.2477 0.6326 500.0068
chain:  1  iteration:  70  parameters:  0.0482 0.2155 0.6698 500.0058
chain:  1  iteration:  80  parameters:  0.0592 0.1971 0.6101 500.054
chain:  1  iteration:  90  parameters:  0.0379 0.2115 0.6569 500.0027
chain:  1  iteration:  100  parameters:  0.0343 0.1978 0.6903 500.0092
chain:  1  iteration:  110  parameters:  0.0386 0.264 0.6412 500.0017
chain:  1  iteration:  120  parameters:  0.0514 0.1758 0.6626 500.0272
chain:  1  iteration:  130  parameters:  0.0558 0.2613 0.5978 500.0098
chain:  1  iteration:  140  parameters:  0.0523 0.2633 0.607 500.0194
chain:  1  iteration:  150  parameters:  0.0586 0.2499 0.5598 500.006
chain:  1  iteration:  160  parameters:  0.0658 0.2747 0.5133 500.0178
chain:  1  iteration:  170  parameters:  0.0536 0.2966 0.5608 500.0001
chain:  1  iteration:  180  parameters:  0.0596 0.3226 0.5326 500.0014
chain:  1  iteration:  190  parameters:  0.0568 0.2861 0.5559 500.0148
chain:  1  iteration:  200  parameters:  0.0496 0.2487 0.6088 500.0012
chain:  2  iteration:  10  parameters:  0.0499 0.2512 0.619 500.0129
chain:  2  iteration:  20  parameters:  0.0515 0.2775 0.5862 500.0071
chain:  2  iteration:  30  parameters:  0.0888 0.2709 0.4641 500.0077
chain:  2  iteration:  40  parameters:  0.0398 0.2471 0.6393 500.008
chain:  2  iteration:  50  parameters:  0.0453 0.1967 0.6801 500.012
chain:  2  iteration:  60  parameters:  0.0479 0.1889 0.6263 500.0109
chain:  2  iteration:  70  parameters:  0.0555 0.1947 0.59 500.0165
chain:  2  iteration:  80  parameters:  0.0556 0.1849 0.616 500.0051
chain:  2  iteration:  90  parameters:  0.0629 0.2343 0.5562 500.0044
chain:  2  iteration:  100  parameters:  0.0446 0.2957 0.6098 500.0145
chain:  2  iteration:  110  parameters:  0.0496 0.2282 0.6275 500.009
chain:  2  iteration:  120  parameters:  0.0419 0.23 0.628 500.003
chain:  2  iteration:  130  parameters:  0.0411 0.212 0.6567 500.0022
chain:  2  iteration:  140  parameters:  0.0421 0.2509 0.6064 500.0081
chain:  2  iteration:  150  parameters:  0.0423 0.2523 0.6228 500.0024
chain:  2  iteration:  160  parameters:  0.0473 0.2189 0.651 500.005
chain:  2  iteration:  170  parameters:  0.0345 0.2256 0.6544 500.0052
chain:  2  iteration:  180  parameters:  0.041 0.1765 0.6751 500.019
chain:  2  iteration:  190  parameters:  0.0446 0.1543 0.6967 500.0161
chain:  2  iteration:  200  parameters:  0.0451 0.2114 0.6256 500.0235
> addPriorConditions <- function(psi){psi[2] + psi[3] < 1}
> MCMC <- bayesGARCH(y, lambda = 100, delta = 500,
+                    control = list(n.chain = 2, l.chain = 200,
+                                   addPriorConditions = addPriorConditions))
chain:  1  iteration:  10  parameters:  0.0418 0.1647 0.6938 500.017
chain:  1  iteration:  20  parameters:  0.0373 0.2278 0.6769 500.0137
chain:  1  iteration:  30  parameters:  0.0403 0.1841 0.7026 500.0095
chain:  1  iteration:  40  parameters:  0.0341 0.1974 0.7136 500.0191
chain:  1  iteration:  50  parameters:  0.0437 0.1844 0.6881 500.0034
chain:  1  iteration:  60  parameters:  0.0296 0.2519 0.693 500.0129
chain:  1  iteration:  70  parameters:  0.0419 0.17 0.7143 500.0066
chain:  1  iteration:  80  parameters:  0.0306 0.2095 0.7261 500.0046
chain:  1  iteration:  90  parameters:  0.0303 0.2363 0.6814 500.0001
chain:  1  iteration:  100  parameters:  0.0402 0.1972 0.6739 500.0306
chain:  1  iteration:  110  parameters:  0.0324 0.2223 0.6892 500.0017
chain:  1  iteration:  120  parameters:  0.0375 0.2387 0.657 500.0184
chain:  1  iteration:  130  parameters:  0.0386 0.2122 0.6721 500.0128
chain:  1  iteration:  140  parameters:  0.051 0.2148 0.6066 500.0318
chain:  1  iteration:  150  parameters:  0.0646 0.2342 0.5599 500.004
chain:  1  iteration:  160  parameters:  0.0513 0.2411 0.5915 500.0097
chain:  1  iteration:  170  parameters:  0.0546 0.3589 0.5126 500.0011
chain:  1  iteration:  180  parameters:  0.0901 0.3275 0.4342 500.0031
chain:  1  iteration:  190  parameters:  0.0533 0.2813 0.5867 500.0005
chain:  1  iteration:  200  parameters:  0.0595 0.2697 0.5207 500.0066
chain:  2  iteration:  10  parameters:  0.0498 0.2439 0.6071 500.0366
chain:  2  iteration:  20  parameters:  0.0631 0.2553 0.5587 500.0086
chain:  2  iteration:  30  parameters:  0.0436 0.1827 0.689 500.0087
chain:  2  iteration:  40  parameters:  0.0412 0.1606 0.704 500.0051
chain:  2  iteration:  50  parameters:  0.0409 0.2078 0.663 500.0177
chain:  2  iteration:  60  parameters:  0.059 0.1848 0.6412 500.0098
chain:  2  iteration:  70  parameters:  0.0622 0.23 0.5504 500.0005
chain:  2  iteration:  80  parameters:  0.0583 0.1872 0.625 500.0063
chain:  2  iteration:  90  parameters:  0.0367 0.3141 0.5929 500.008
chain:  2  iteration:  100  parameters:  0.0576 0.2113 0.633 500.0277
chain:  2  iteration:  110  parameters:  0.0505 0.2083 0.6078 500
chain:  2  iteration:  120  parameters:  0.0491 0.1966 0.6778 500.0036
chain:  2  iteration:  130  parameters:  0.0351 0.2449 0.6538 500.0006
chain:  2  iteration:  140  parameters:  0.0483 0.1964 0.634 500.005
chain:  2  iteration:  150  parameters:  0.0398 0.2056 0.6834 500.0157
chain:  2  iteration:  160  parameters:  0.0468 0.1898 0.6829 500.0036
chain:  2  iteration:  170  parameters:  0.0372 0.2515 0.6407 500.0018
chain:  2  iteration:  180  parameters:  0.0312 0.1982 0.7218 500.0054
chain:  2  iteration:  190  parameters:  0.0367 0.1801 0.7112 500.0004
chain:  2  iteration:  200  parameters:  0.0308 0.2253 0.7052 500.004
> install.packages("fGarch")
Installing package into ‘C:/Users/ADMIN/Documents/R/win-library/3.5’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.5/fGarch_3042.83.1.zip'
Content type 'application/zip' length 607349 bytes (593 KB)
downloaded 593 KB

package ‘fGarch’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
C:\Users\ADMIN\AppData\Local\Temp\Rtmp6VdpoA\downloaded_packages
> library(fGarch)

Attaching package: ‘fGarch’

The following object is masked _by_ ‘.GlobalEnv’:

    dem2gbp

Warning message:
package ‘fGarch’ was built under R version 3.5.3
> x = as.timeSeries(data(LPP2005REC))
> fit = garchFit(LPP40 ~ garch(1, 1), data = 100*x, trace = FALSE)
> volatility = volatility(fit, type = "sigma")
> head(volatility)
[1] 0.2805121 0.2674818 0.2608402
[4] 0.2645372 0.2603407 0.2558201
> class(volatility)
[1] "numeric"
> volatility = volatility(fit, type = "h")
> head(volatility)
[1] 0.07868705 0.07154651 0.06803759
[4] 0.06997992 0.06777728 0.06544394
> class(volatility)
[1] "numeric"
> volatility = slot(fit, "sigma.t")
> head(volatility)
[1] 0.2805121 0.2674818 0.2608402
[4] 0.2645372 0.2603407 0.2558201
> class(volatility)
[1] "numeric"
> volatility = slot(fit, "sigma.t")
> head(volatility)
[1] 0.2805121 0.2674818 0.2608402
[4] 0.2645372 0.2603407 0.2558201
> class(volatility)
[1] "numeric"
> volatility = slot(fit, "h.t")
> head(volatility)
[1] 0.07868705 0.07154651 0.06803759
[4] 0.06997992 0.06777728 0.06544394
> par(mfrow = c(2, 2))
> set.seed(1953)
> r = rsged(n = 1000)
> plot(r, type = "l", main = "sged", col = "steelblue")
> hist(r, n = 25, probability = TRUE, border = "white", col = "steelblue")
> box()
> x = seq(min(r), max(r), length = 201)
> lines(x, dsged(x), lwd = 2)
> plot(sort(r), (1:1000/1000), main = "Probability", col = "steelblue",
+      ylab = "Probability")
> lines(x, psged(x), lwd = 2)
> round(qsged(psged(q = seq(-1, 5, by = 1))), digits = 6)
[1] -1  0  1  2  3  4  5
> plot(sort(r), (1:1000/1000), main = "Probability", col = "steelblue",
+      ylab = "Probability")
> lines(x, psged(x), lwd = 2)
> round(qsged(psged(q = seq(-1, 5, by = 1))), digits = 6)
[1] -1  0  1  2  3  4  5
> set.seed(123)
> fit = garchFit(~ garch(1, 1), data = garchSim(), trace = FALSE)
> fit

Title:
 GARCH Modelling

Call:
 garchFit(formula = ~garch(1, 1), data = garchSim(), trace = FALSE)

Mean and Variance Equation:
 data ~ garch(1, 1)
<environment: 0x0000000012ba1338>
 [data = garchSim()]

Conditional Distribution:
 norm

Coefficient(s):
         mu        omega       alpha1
-1.5658e-05   3.1101e-06   2.8879e-01
      beta1
 4.0817e-01

Std. Errors:
 based on Hessian

Error Analysis:
         Estimate  Std. Error  t value
mu     -1.566e-05   2.637e-04   -0.059
omega   3.110e-06   1.874e-06    1.660
alpha1  2.888e-01   1.808e-01    1.597
beta1   4.082e-01   2.777e-01    1.470
       Pr(>|t|)
mu        0.953
omega     0.097 .
alpha1    0.110
beta1     0.142
---
Signif. codes:
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’
  0.1 ‘ ’ 1

Log Likelihood:
 440.3751    normalized:  4.403751

Description:
 Sat Apr 20 11:26:49 2019 by user: ADMIN

> predict(fit, n.ahead = 10)
    meanForecast   meanError
1  -1.565789e-05 0.004533623
2  -1.565789e-05 0.004175561
3  -1.565789e-05 0.003906646
4  -1.565789e-05 0.003707707
5  -1.565789e-05 0.003562490
6  -1.565789e-05 0.003457674
7  -1.565789e-05 0.003382702
8  -1.565789e-05 0.003329451
9  -1.565789e-05 0.003291827
10 -1.565789e-05 0.003265349
   standardDeviation
1        0.004533623
2        0.004175561
3        0.003906646
4        0.003707707
5        0.003562490
6        0.003457674
7        0.003382702
8        0.003329451
9        0.003291827
10       0.003265349
> predict(fit, n.ahead = 10,mse="uncond")
    meanForecast  meanError
1  -1.565789e-05 0.00305647
2  -1.565789e-05 0.00305647
3  -1.565789e-05 0.00305647
4  -1.565789e-05 0.00305647
5  -1.565789e-05 0.00305647
6  -1.565789e-05 0.00305647
7  -1.565789e-05 0.00305647
8  -1.565789e-05 0.00305647
9  -1.565789e-05 0.00305647
10 -1.565789e-05 0.00305647
   standardDeviation
1        0.004533623
2        0.004175561
3        0.003906646
4        0.003707707
5        0.003562490
6        0.003457674
7        0.003382702
8        0.003329451
9        0.003291827
10       0.003265349
> predict(fit, n.ahead = 10, plot=TRUE, crit_val=2)
    meanForecast   meanError
1  -1.565789e-05 0.004533623
2  -1.565789e-05 0.004175561
3  -1.565789e-05 0.003906646
4  -1.565789e-05 0.003707707
5  -1.565789e-05 0.003562490
6  -1.565789e-05 0.003457674
7  -1.565789e-05 0.003382702
8  -1.565789e-05 0.003329451
9  -1.565789e-05 0.003291827
10 -1.565789e-05 0.003265349
   standardDeviation lowerInterval
1        0.004533623  -0.009082904
2        0.004175561  -0.008366780
3        0.003906646  -0.007828950
4        0.003707707  -0.007431071
5        0.003562490  -0.007140637
6        0.003457674  -0.006931006
7        0.003382702  -0.006781061
8        0.003329451  -0.006674559
9        0.003291827  -0.006599312
10       0.003265349  -0.006546355
   upperInterval
1    0.009051588
2    0.008335464
3    0.007797635
4    0.007399756
5    0.007109322
6    0.006899691
7    0.006749746
8    0.006643243
9    0.006567996
10   0.006515039
> set.seed(321)
> fit2 = garchFit(~ garch(1, 1), data = garchSim(), trace = FALSE, cond.dist="sged")
https://www.mathclasstutor.online

https://www.mathclasstutor.online

> predict(fit2,n.ahead=20,plot=TRUE)
    meanForecast   meanError
1  -0.0001011762 0.002792102
2  -0.0001011762 0.002697756
3  -0.0001011762 0.002657046
4  -0.0001011762 0.002639724
5  -0.0001011762 0.002632400
6  -0.0001011762 0.002629311
7  -0.0001011762 0.002628010
8  -0.0001011762 0.002627462
9  -0.0001011762 0.002627231
10 -0.0001011762 0.002627134
11 -0.0001011762 0.002627093
12 -0.0001011762 0.002627076
13 -0.0001011762 0.002627069
14 -0.0001011762 0.002627066
15 -0.0001011762 0.002627064
16 -0.0001011762 0.002627064
17 -0.0001011762 0.002627064
18 -0.0001011762 0.002627064
19 -0.0001011762 0.002627064
20 -0.0001011762 0.002627064

   StandardDeviation lowerInterval

1        0.002792102  -0.005672090
2        0.002697756  -0.005483848
3        0.002657046  -0.005402622
4        0.002639724  -0.005368060
5        0.002632400  -0.005353446
6        0.002629311  -0.005347283
7        0.002628010  -0.005344686
8        0.002627462  -0.005343593
9        0.002627231  -0.005343133
10       0.002627134  -0.005342940
11       0.002627093  -0.005342858
12       0.002627076  -0.005342824
13       0.002627069  -0.005342809
14       0.002627066  -0.005342803
15       0.002627064  -0.005342801
16       0.002627064  -0.005342800
17       0.002627064  -0.005342799
18       0.002627064  -0.005342799
19       0.002627064  -0.005342799
20       0.002627064  -0.005342799

   upperInterval

1    0.005183724
2    0.005005147
3    0.004928091
4    0.004895303
5    0.004881439
6    0.004875592
7    0.004873130
8    0.004872093
9    0.004871656
10   0.004871472
11   0.004871395
12   0.004871363
13   0.004871349
14   0.004871343
15   0.004871341
16   0.004871340
17   0.004871339
18   0.004871339
19   0.004871339
20   0.004871339

Friday, 19 April 2019

How to block maxima for Stock losses

 How to block maxima for Stock losses

R is free software
> library ( evir )
Warning messages:
1: package ‘fOptions’ 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.2
5: package ‘evir’ was built under R version 3.5.2
> data ( siemens )
> SieLoss <− −100.0 ∗ siemens
> SieGEV <− gev ( SieLoss , block = " semester " )
Error in gev(SieLoss, block = " semester ") : unknown time period
> library(evir)
> SieGEV <− gev ( SieLoss , block = " semester " )
Error in gev(SieLoss, block = " semester ") : unknown time period
> out <- gev(bmw, "month")
Error in gev(bmw, "month") : object 'bmw' not found
> require(graphics)
> fr <- function(x) {   ## Rosenbrock Banana function
+     x1 <- x[1]
+     x2 <- x[2]
+     100 * (x2 - x1 * x1)^2 + (1 - x1)^2
+ }
> grr <- function(x) { ## Gradient of 'fr'
+     x1 <- x[1]
+     x2 <- x[2]
+     c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1),
+       200 *      (x2 - x1 * x1))
+ }
> optim(c(-1.2,1), fr)
$`par`
[1] 1.000260 1.000506

$value
[1] 8.825241e-08

$counts
function gradient
     195       NA

$convergence
[1] 0

$message
NULL

> (res <- optim(c(-1.2,1), fr, grr, method = "BFGS"))
$`par`
[1] 1 1

$value
[1] 9.594956e-18

$counts
function gradient
     110       43

$convergence
[1] 0

$message
NULL

> optimHess(res$par, fr, grr)
          [,1] [,2]
[1,]  802.0004 -400
[2,] -400.0000  200
> optim(c(-1.2,1), fr, NULL, method = "BFGS", hessian = TRUE)
$`par`
[1] 0.9998044 0.9996084

$value
[1] 3.827383e-08

$counts
function gradient
     118       38

$convergence
[1] 0

$message
NULL

$hessian
          [,1]      [,2]
[1,]  801.6881 -399.9218
[2,] -399.9218  200.0000

> optim(c(-1.2,1), fr, grr, method = "CG")
$`par`
[1] -0.7648373  0.5927588

$value
[1] 3.106579

$counts
function gradient
     402      101

$convergence
[1] 1

$message
NULL

> optim(c(-1.2,1), fr, grr, method = "CG", control = list(type = 2))
$`par`
[1] 0.9944093 0.9888229

$value
[1] 3.123777e-05

$counts
function gradient
     385      101

$convergence
[1] 1

$message
NULL

> optim(c(-1.2,1), fr, grr, method = "L-BFGS-B")
$`par`
[1] 0.9999997 0.9999995

$value
[1] 2.267577e-13

$counts
function gradient
      47       47

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

> flb <- function(x)
+ { p <- length(x); sum(c(1, rep(4, p-1)) * (x - c(1, x[-p])^2)^2) }
> optim(rep(3, 25), flb, NULL, method = "L-BFGS-B",
+       lower = rep(2, 25), upper = rep(4, 25))
$`par`
 [1] 2.000000 2.000000 2.000000 2.000000
 [5] 2.000000 2.000000 2.000000 2.000000
 [9] 2.000000 2.000000 2.000000 2.000000
[13] 2.000000 2.000000 2.000000 2.000000
[17] 2.000000 2.000000 2.000000 2.000000
[21] 2.000000 2.000000 2.000000 2.109093
[25] 4.000000

$value
[1] 368.1059

$counts
function gradient
       6        6

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

>
> fw <- function (x)
+     10*sin(0.3*x)*sin(1.3*x^2) + 0.00001*x^4 + 0.2*x+80
>
> plot(fw, -50, 50, n = 1000, main = "optim() minimising 'wild function'")
>
> res <- optim(50, fw, method = "SANN",
+              control = list(maxit = 20000, temp = 20, parscale = 20))
> res
$`par`
[1] -15.8146

$value
[1] 67.4703

$counts
function gradient
   20000       NA

$convergence
[1] 0

$message
NULL

> (r2 <- optim(res$par, fw, method = "BFGS"))
$`par`
[1] -15.81515

$value
[1] 67.46773

$counts
function gradient
      16        3

$convergence
[1] 0

$message
NULL

> points(r2$par, r2$value, pch = 8, col = "red", cex = 2)
> library(stats)
> eurodistmat <- as.matrix(eurodist)
>
> distance <- function(sq) {  # Target function
+     sq2 <- embed(sq, 2)
+     sum(eurodistmat[cbind(sq2[,2], sq2[,1])])
+ }
> genseq <- function(sq) {  # Generate new candidate sequence
+     idx <- seq(2, NROW(eurodistmat)-1)
+     changepoints <- sample(idx, size = 2, replace = FALSE)
+     tmp <- sq[changepoints[1]]
+     sq[changepoints[1]] <- sq[changepoints[2]]
+     sq[changepoints[2]] <- tmp
+     sq
+ }
> sq <- c(1:nrow(eurodistmat), 1)
> loc <- -cmdscale(eurodist, add = TRUE)$points
> x <- loc[,1]; y <- loc[,2]
> s <- seq_len(nrow(eurodistmat))
> tspinit <- loc[sq,]
>
> plot(x, y, type = "n", asp = 1, xlab = "", ylab = "",
+      main = "initial solution of traveling salesman problem", axes = FALSE)
> arrows(tspinit[s,1], tspinit[s,2], tspinit[s+1,1], tspinit[s+1,2],
+        angle = 10, col = "green")
> text(x, y, labels(eurodist), cex = 0.8)
> set.seed(123)
> res <- optim(sq, distance, genseq, method = "SANN",
+              control = list(maxit = 30000, temp = 2000, trace = TRUE,
+                             REPORT = 500))
sann objective function values
initial       value 29625.000000
iter     5000 value 13585.000000
iter    10000 value 13092.000000
iter    15000 value 13063.000000
iter    20000 value 12919.000000
iter    25000 value 12907.000000
iter    29999 value 12842.000000
final         value 12842.000000
sann stopped after 29999 iterations
> tspres <- loc[res$par,]
> plot(x, y, type = "n", asp = 1, xlab = "", ylab = "",
+      main = "optim() 'solving' traveling salesman problem", axes = FALSE)
> arrows(tspres[s,1], tspres[s,2], tspres[s+1,1], tspres[s+1,2],
+        angle = 10, col = "red")
> text(x, y, labels(eurodist), cex = 0.8)
> system.time(rO <- optimize(function(x) (x-pi)^2, c(0, 10)))
   user  system elapsed
      0       0       0
> system.time(ro <- optim(1, function(x) (x-pi)^2, control=list(warn.1d.NelderMead = FALSE)))
   user  system elapsed
      0       0       0
> rO$minimum - pi # 0 (perfect), on one platform
[1] 0
> ro$par - pi     # ~= 1.9e-4    on one platform
[1] -0.0001864036
> utils::str(ro)
List of 5
 $ par        : num 3.14
 $ value      : num 3.47e-08
 $ counts     : Named int [1:2] 32 NA
  ..- attr(*, "names")= chr [1:2] "function" "gradient"
 $ convergence: int 0
 $ message    : NULL
> library("evir", lib.loc="~/R/win-library/3.5")
> SieGEV <− gev ( SieLoss , block = " semester " )
Error in gev(SieLoss, block = " semester ") : unknown time period
> SieGEV <− gev ( SieLoss , block = " month " )
Error in gev(SieLoss, block = " month ") : unknown time period
> data(bmw)
> out <- gev(bmw, "month")
> out
$`n.all`
[1] 6146

$n
[1] 283

$data
  [1] 0.047704097 0.040072550 0.057006818
  [4] 0.042967868 0.022428128 0.040120865
  [7] 0.030498439 0.025262456 0.057779537
 [10] 0.039382379 0.031748698 0.021695014
 [13] 0.085708790 0.017057713 0.073563768
 [16] 0.018876364 0.033584741 0.058143393
 [19] 0.042326924 0.009595910 0.032909034
 [22] 0.051069736 0.043738236 0.050509198
 [25] 0.066278121 0.040930359 0.036095789
 [28] 0.029984390 0.020905075 0.032823886
 [31] 0.026219918 0.047408878 0.047743616
 [34] 0.018994386 0.027733998 0.031215494
 [37] 0.027928237 0.011354940 0.040061610
 [40] 0.015156798 0.016943433 0.015530942
 [43] 0.022243069 0.012447720 0.010386600
 [46] 0.029017005 0.015888185 0.012732589
 [49] 0.010195424 0.018930242 0.027147614
 [52] 0.017685602 0.011002997 0.026742729
 [55] 0.023126669 0.010324201 0.016317671
 [58] 0.022688656 0.018503749 0.029676230
 [61] 0.017803681 0.028395043 0.011157337
 [64] 0.020159906 0.008803704 0.027067868
 [67] 0.013293438 0.020227519 0.028373642
 [70] 0.017675785 0.017558952 0.015440693
 [73] 0.020048081 0.010208282 0.038070120
 [76] 0.029718365 0.015372579 0.010795260
 [79] 0.008982199 0.030438497 0.068411993
 [82] 0.010744914 0.033322322 0.011913702
 [85] 0.013162712 0.027771744 0.019664256
 [88] 0.040528019 0.018800955 0.033664422
 [91] 0.024044674 0.027267717 0.016999328
 [94] 0.069930216 0.025860543 0.022775877
 [97] 0.035147558 0.029596241 0.023585409
[100] 0.012656443 0.020927520 0.040695897
[103] 0.019978081 0.031988439 0.025409997
[106] 0.021238736 0.045243354 0.039014446
[109] 0.013019157 0.044288785 0.027986842
[112] 0.028975106 0.015568407 0.014595571
[115] 0.017141573 0.014628069 0.028136965
[118] 0.022023821 0.022388485 0.017276267
[121] 0.018211336 0.028117197 0.029418132
[124] 0.017657790 0.038427474 0.032953268
[127] 0.015801683 0.037267309 0.034521814
[130] 0.030317974 0.015823230 0.040316266
[133] 0.021733491 0.020395699 0.021325364
[136] 0.018883156 0.021034092 0.016287900
[139] 0.026655237 0.025479085 0.016252641
[142] 0.030026651 0.029463070 0.054088773
[145] 0.030317974 0.050467667 0.090348294
[148] 0.028732294 0.016799955 0.011537384
[151] 0.029699096 0.066342789 0.019648619
[154] 0.051285118 0.041242959 0.033119684
[157] 0.088410957 0.030153038 0.034251807
[160] 0.011888294 0.055775683 0.046194152
[163] 0.028076501 0.034191365 0.044215223
[166] 0.066280459 0.044012805 0.027359104
[169] 0.035537302 0.043172876 0.117191789
[172] 0.053856606 0.043313442 0.046766994
[175] 0.029291783 0.033416901 0.036793770
[178] 0.029772202 0.016362673 0.071609532
[181] 0.085028505 0.028179646 0.013151022
[184] 0.016368476 0.018105128 0.024107501
[187] 0.040991244 0.026881436 0.024632454
[190] 0.012288641 0.024679312 0.028208721
[193] 0.018634561 0.021706789 0.024693705
[196] 0.057176834 0.015042360 0.012608702
[199] 0.017940725 0.029723723 0.036431458
[202] 0.026984661 0.017313206 0.057152069
[205] 0.051660609 0.041497805 0.033114946
[208] 0.026184444 0.013923817 0.026675428
[211] 0.029396614 0.019142057 0.021650150
[214] 0.059098050 0.039794958 0.067674756
[217] 0.088632045 0.062040653 0.041325860
[220] 0.016045377 0.042388410 0.040997961
[223] 0.018651503 0.018750425 0.019290632
[226] 0.044707716 0.005662407 0.009850060
[229] 0.023593976 0.014741345 0.018630817
[232] 0.018749337 0.016068787 0.018174172
[235] 0.035958313 0.007357140 0.024359967
[238] 0.025035790 0.055598597 0.035407724
[241] 0.020823390 0.011553667 0.017623795
[244] 0.023831384 0.020582821 0.029107995
[247] 0.012860773 0.021955587 0.068605893
[250] 0.016159641 0.014116384 0.037276749
[253] 0.079474594 0.066689116 0.016613801
[256] 0.017747906 0.024150295 0.018389544
[259] 0.016504432 0.019881371 0.035009410
[262] 0.018259789 0.009618542 0.030405515
[265] 0.013192803 0.030913346 0.021380795
[268] 0.011232361 0.019076147 0.022407876
[271] 0.031367270 0.025284978 0.024784416
[274] 0.016912398 0.009060769 0.022973983
[277] 0.033093688 0.015351852 0.016209831
[280] 0.026569373 0.017001017 0.015623468
[283] 0.019133955

$block
[1] "month"

$par.ests
        xi      sigma         mu
0.22084471 0.01016219 0.02104994

$par.ses
          xi        sigma           mu
0.0538218570 0.0005249399 0.0006934577

$varcov
              [,1]          [,2]          [,3]
[1,]  2.896792e-03 -2.617132e-06 -1.249425e-05
[2,] -2.617132e-06  2.755619e-07  2.157608e-07
[3,] -1.249425e-05  2.157608e-07  4.808836e-07

$converged
[1] 0

$nllh.final
[1] -816.0357

attr(,"class")
[1] "gev"
>
>
>
> out<- gev(bmw,"month")
> out <- gev(bmw, 100)
> SieGEV <− gev( SieLoss,block ="month")
> plot (SieGEV$data,type = " h " , col = " blue ", xlab = " " ,
+         ylab = " Block Maxima" ,
+         main = "Maximum BiannualLossesof Siemens " )
Error in plot.xy(xy, type, ...) : invalid plot type ' '
In addition: Warning message:
In plot.xy(xy, type, ...) :
  plot type ' h ' will be truncated to first character
> plot (SieGEV$data,type = " h ", col = " blue ", xlab = " " ,ylab = " Block Maxima" , main = "Maximum BiannualLossesof Siemens ")
Error in plot.xy(xy, type, ...) : invalid plot type ' '
In addition: Warning message:
In plot.xy(xy, type, ...) :
  plot type ' h ' will be truncated to first character
> plot (SieGEV$data,type = "h", col = "blue", xlab = " ,ylab = " Block Maxima" , main = "Maximum Biannual Losses of Siemens")
Error: unexpected symbol in "plot (SieGEV$data,type = "h", col = "blue", xlab = " ,ylab = " Block"
> SieGEV
$`n.all`
[1] 6146

$n
[1] 283

$data
  [1]  1.6536690  3.4026271  3.5896764
  [4]  2.4363154  2.2190569  2.6621690
  [7]  3.0927687  3.5197198  2.4184353
 [10]  2.3467744  2.1210311  1.5029506
 [13]  2.1874949  3.4257494  2.0331069
 [16]  2.8475728  2.2709367  2.0191972
 [19]  1.4117882  3.1748698  1.6491702
 [22]  2.2267549  3.1872617  1.9841921
 [25]  1.5604298  1.4921810  0.9151478
 [28]  0.5331159  1.3761685  1.8930788
 [31]  2.7684999  0.9991000  1.0722813
 [34]  1.8397365  1.5435808  1.0640662
 [37]  0.9368628  1.0730908  1.9506175
 [40]  1.3471241  3.5896546  1.6234766
 [43]  1.7286863  1.2846423  0.9982115
 [46]  2.1193397  1.0924478  2.4508885
 [49]  2.9427001  1.6106459  0.8926898
 [52]  1.2969331  1.2121361  0.5624068
 [55]  2.1739987  1.6633335  0.7278987
 [58]  0.9950331  0.5942293  0.7528266
 [61]  0.9508788  0.8940779  1.6221919
 [64]  1.0400094  2.3063940  1.6047642
 [67]  0.5065439  0.6480379  1.1514918
 [70]  1.3880849  0.6595562  1.8634080
 [73]  1.0929071  1.7669791  1.7094433
 [76]  0.8688152  2.7912846  2.6773598
 [79]  0.8360885  1.4056001  0.9356272
 [82]  1.1575692  0.5542542  2.0636656
 [85]  1.1870880  1.8915038  1.6469543
 [88]  1.6116384  4.9999547  1.1622840
 [91]  0.7788521  0.6306682  0.6920443
 [94]  1.5674302  0.8240188  1.4173466
 [97]  1.2407107  1.7079834  2.2397353
[100]  3.3175862  3.3228516  0.5922183
[103]  1.7978998  0.7527155  1.3144137
[106]  3.5107748  2.7424887  3.9298765
[109]  1.7241806  0.6972140  1.0398707
[112]  1.1096998  3.5335848  1.7425417
[115]  1.0779134  1.6385415  1.1560822
[118]  2.2122571  4.4886451  2.5135973
[121]  1.7623495  0.7214460  3.1989716
[124]  0.9447139  1.4484497  1.4198687
[127]  3.2918569  1.8343128  1.9322879
[130]  2.3835173  1.1528403  0.7318097
[133]  1.5378692  3.2633218  0.8741717
[136]  0.8289248  1.5479779  1.1374279
[139]  1.7396470  1.4001306  2.0545854
[142]  1.2002326  0.9077941  1.7542311
[145]  1.5805540  1.2834196  2.8254668
[148]  2.8259337  3.0883472  1.0949014
[151]  1.6965534  2.4391453  3.4268457
[154]  1.4953550  2.4938948  2.6579638
[157]  3.5718083  4.7905600  1.1379924
[160]  1.3513719  6.6021101  4.0637646
[163]  3.8652154  2.5001302  2.8358865
[166]  1.4556298  2.6126305  1.7518696
[169]  6.3867284  3.1010237  8.6568016
[172]  3.7944932  2.8393075  3.4635497
[175]  1.4144507  2.4259950  5.6995481
[178]  2.6472134  1.9652938  7.2320662
[181]  4.7829088  1.1770863  1.7784756
[184]  1.0294897  8.1770542  3.3198069
[187]  3.2595478  1.3297028  2.5196294
[190]  2.0831702  1.0382349  2.0367303
[193]  3.0934190  1.9465795  1.3685825
[196]  1.4377849  2.3776127  0.7466900
[199]  1.5325970  1.0976676  0.9890283
[202]  1.8210658  2.2436317 12.0111624
[205]  2.6916164  2.5308099  2.7487583
[208]  2.8735689  1.5320721  2.1862414
[211]  1.9762419  1.5716073  1.2565746
[214]  6.0060331  2.9175489  2.8963547
[217]  3.4198861  2.3000384  1.3762774
[220]  1.4864023  2.7354883  1.7335580
[223]  0.9072227  2.3636616  1.2646683
[226]  9.3806420  1.0966208  0.9709587
[229]  0.9623170  0.7572981  1.2196383
[232]  0.9762859  2.0181874  1.1351194
[235]  0.6193797  1.3728856  2.7488846
[238]  1.8836173  1.4787700  3.3931668
[241]  2.2802854  1.0756865  2.2783174
[244]  1.2495070  2.8341671  1.2006069
[247]  1.1539554  1.1361914  1.8208750
[250]  2.1920309  1.9187432  0.9584072
[253]  4.9052494  2.5863511  2.2497957
[256]  2.0670434  2.0282303  2.0100492
[259]  2.0702936  3.6473998  1.3965681
[262]  1.5479044  3.4507620  3.2426724
[265]  1.4998751  1.1565900  1.2561791
[268]  0.8417013  2.0926198  2.0305266
[271]  1.6009918  1.1278315  1.9391189
[274]  0.9086656  3.3403239  2.9323615
[277]  0.5704793  1.4510533  1.7931181
[280]  1.1174629  1.5387715  1.4401030
[283]  2.1357596

$block
[1] "month"

$par.ests
       xi     sigma        mu
0.2591183 0.7181868 1.4309933

$par.ses
        xi      sigma         mu
0.05438777 0.04055958 0.04922595

$varcov
              [,1]         [,2]          [,3]
[1,]  0.0029580294 -0.000196826 -0.0008900357
[2,] -0.0001968260  0.001645080  0.0012402593
[3,] -0.0008900357  0.001240259  0.0024231940

$converged
[1] 0

$nllh.final
[1] 394.6513

attr(,"class")
[1] "gev"
> plot (SieGEV$data,type = "h", col = "blue", xlab = " ,ylab = " Block Maxima" , main = "Maximum Biannual Losses of Siemens")
Error: unexpected symbol in "plot (SieGEV$data,type = "h", col = "blue", xlab = " ,ylab = " Block"
> plot ( SieGEV$data , type = " h " , col = " blue " , xlab = " " ,
+         ylab = " Block Maxima" ,
+         main = "Maximum Biannual Losses of Siemens " )
Error in plot.xy(xy, type, ...) : invalid plot type ' '
In addition: Warning message:
In plot.xy(xy, type, ...) :
  plot type ' h ' will be truncated to first character
> plot(out)

Make a plot selection (or 0 to exit):

1: plot: Scatterplot of Residuals
2: plot: QQplot of Residuals

Selection:
Enter an item from the menu, or 0 to exit
Selection:
Enter an item from the menu, or 0 to exit
Selection: 1

Make a plot selection (or 0 to exit):

1: plot: Scatterplot of Residuals
2: plot: QQplot of Residuals

Selection: 2

Make a plot selection (or 0 to exit):

1: plot: Scatterplot of Residuals
2: plot: QQplot of Residuals

Selection: 3
Enter an item from the menu, or 0 to exit
Selection: 0
> data(danish)
> qplot(danish)
> plot(out)

Make a plot selection (or 0 to exit):

1: plot: Scatterplot of Residuals
2: plot: QQplot of Residuals

Selection: 1

Make a plot selection (or 0 to exit):

1: plot: Scatterplot of Residuals
2: plot: QQplot of Residuals

Selection: 2

Make a plot selection (or 0 to exit):

1: plot: Scatterplot of Residuals
2: plot: QQplot of Residuals

Selection: 3
Enter an item from the menu, or 0 to exit
Selection: 4
Enter an item from the menu, or 0 to exit
Selection: 0
> data(bmw)
> out <- pot(-bmw, ne = 200)
> decluster(out$data, 30)
Declustering picture...
Data reduced from 200 to 70
 [1] 0.05525992 0.05084288 0.06689143
 [4] 0.06887800 0.04688708 0.10617520
 [7] 0.05865478 0.04071008 0.02896766
[10] 0.02913723 0.04122546 0.03094040
[13] 0.02659780 0.03465775 0.02606819
[16] 0.02710781 0.06955422 0.03859640
[19] 0.03845126 0.06807869 0.02732490
[22] 0.04194781 0.02563799 0.04957226
[25] 0.03744841 0.02652559 0.03021099
[28] 0.04541566 0.02695412 0.03051167
[31] 0.03984488 0.02648610 0.04525858
[34] 0.03938106 0.03236528 0.05561317
[37] 0.03191673 0.04609111 0.06519183
[40] 0.03668970 0.03540193 0.04491698
[43] 0.04880244 0.04565154 0.04148352
[46] 0.02619353 0.10852160 0.02954286
[49] 0.03225168 0.14061565 0.02895132
[52] 0.02834847 0.03360120 0.07529138
[55] 0.04328683 0.03405556 0.02677080
[58] 0.03221330 0.10577463 0.02764203
[61] 0.05806183 0.03004121 0.03973532
[64] 0.02892858 0.02688352 0.04001846
[67] 0.03412405 0.03232889 0.03123302
[70] 0.03490287
attr(,"times")
 [1] "1973-02-05 05:30:00 IST"
 [2] "1973-07-26 05:30:00 IST"
 [3] "1973-11-14 05:30:00 IST"
 [4] "1973-12-31 05:30:00 IST"
 [5] "1974-05-09 05:30:00 IST"
 [6] "1974-07-05 05:30:00 IST"
 [7] "1974-09-30 05:30:00 IST"
 [8] "1975-05-27 05:30:00 IST"
 [9] "1976-04-06 05:30:00 IST"
[10] "1976-07-28 05:30:00 IST"
[11] "1976-10-12 05:30:00 IST"
[12] "1976-11-15 05:30:00 IST"
[13] "1977-02-18 05:30:00 IST"
[14] "1977-07-08 05:30:00 IST"
[15] "1977-12-06 05:30:00 IST"
[16] "1978-06-21 05:30:00 IST"
[17] "1979-07-11 05:30:00 IST"
[18] "1980-01-21 05:30:00 IST"
[19] "1980-03-18 05:30:00 IST"
[20] "1980-06-20 05:30:00 IST"
[21] "1980-09-08 05:30:00 IST"
[22] "1981-01-27 05:30:00 IST"
[23] "1981-04-29 05:30:00 IST"
[24] "1981-06-26 05:30:00 IST"
[25] "1981-08-25 05:30:00 IST"
[26] "1982-04-27 05:30:00 IST"
[27] "1982-08-17 05:30:00 IST"
[28] "1982-09-27 05:30:00 IST"
[29] "1983-01-24 05:30:00 IST"
[30] "1983-07-12 05:30:00 IST"
[31] "1984-02-09 05:30:00 IST"
[32] "1984-03-12 05:30:00 IST"
[33] "1984-07-05 05:30:00 IST"
[34] "1984-10-19 05:30:00 IST"
[35] "1985-03-25 05:30:00 IST"
[36] "1985-07-10 05:30:00 IST"
[37] "1985-09-13 05:30:00 IST"
[38] "1985-11-05 05:30:00 IST"
[39] "1986-02-27 05:30:00 IST"
[40] "1986-04-28 05:30:00 IST"
[41] "1986-06-03 05:30:00 IST"
[42] "1986-07-04 05:30:00 IST"
[43] "1986-09-05 05:30:00 IST"
[44] "1987-02-04 05:30:00 IST"
[45] "1987-04-27 05:30:00 IST"
[46] "1987-06-02 05:30:00 IST"
[47] "1987-11-09 05:30:00 IST"
[48] "1988-03-25 05:30:00 IST"
[49] "1988-05-11 05:30:00 IST"
[50] "1989-10-16 05:30:00 IST"
[51] "1990-01-04 05:30:00 IST"
[52] "1990-02-26 05:30:00 IST"
[53] "1990-04-23 05:30:00 IST"
[54] "1990-09-25 05:30:00 IST"
[55] "1991-01-14 05:30:00 IST"
[56] "1991-02-26 05:30:00 IST"
[57] "1991-05-17 05:30:00 IST"
[58] "1991-06-28 05:30:00 IST"
[59] "1991-08-19 05:30:00 IST"
[60] "1992-05-13 05:30:00 IST"
[61] "1992-09-24 05:30:00 IST"
[62] "1993-01-25 05:30:00 IST"
[63] "1993-05-14 05:30:00 IST"
[64] "1993-07-27 05:30:00 IST"
[65] "1993-11-22 05:30:00 IST"
[66] "1994-06-20 05:30:00 IST"
[67] "1994-11-23 05:30:00 IST"
[68] "1995-03-09 05:30:00 IST"
[69] "1995-09-22 05:30:00 IST"
[70] "1995-10-23 05:30:00 IST"
> out <- pot(-bmw, ne = 200)
> decluster(out$data, 30)
Declustering picture...
Data reduced from 200 to 70
 [1] 0.05525992 0.05084288 0.06689143
 [4] 0.06887800 0.04688708 0.10617520
 [7] 0.05865478 0.04071008 0.02896766
[10] 0.02913723 0.04122546 0.03094040
[13] 0.02659780 0.03465775 0.02606819
[16] 0.02710781 0.06955422 0.03859640
[19] 0.03845126 0.06807869 0.02732490
[22] 0.04194781 0.02563799 0.04957226
[25] 0.03744841 0.02652559 0.03021099
[28] 0.04541566 0.02695412 0.03051167
[31] 0.03984488 0.02648610 0.04525858
[34] 0.03938106 0.03236528 0.05561317
[37] 0.03191673 0.04609111 0.06519183
[40] 0.03668970 0.03540193 0.04491698
[43] 0.04880244 0.04565154 0.04148352
[46] 0.02619353 0.10852160 0.02954286
[49] 0.03225168 0.14061565 0.02895132
[52] 0.02834847 0.03360120 0.07529138
[55] 0.04328683 0.03405556 0.02677080
[58] 0.03221330 0.10577463 0.02764203
[61] 0.05806183 0.03004121 0.03973532
[64] 0.02892858 0.02688352 0.04001846
[67] 0.03412405 0.03232889 0.03123302
[70] 0.03490287
attr(,"times")
 [1] "1973-02-05 05:30:00 IST"
 [2] "1973-07-26 05:30:00 IST"
 [3] "1973-11-14 05:30:00 IST"
 [4] "1973-12-31 05:30:00 IST"
 [5] "1974-05-09 05:30:00 IST"
 [6] "1974-07-05 05:30:00 IST"
 [7] "1974-09-30 05:30:00 IST"
 [8] "1975-05-27 05:30:00 IST"
 [9] "1976-04-06 05:30:00 IST"
[10] "1976-07-28 05:30:00 IST"
[11] "1976-10-12 05:30:00 IST"
[12] "1976-11-15 05:30:00 IST"
[13] "1977-02-18 05:30:00 IST"
[14] "1977-07-08 05:30:00 IST"
[15] "1977-12-06 05:30:00 IST"
[16] "1978-06-21 05:30:00 IST"
[17] "1979-07-11 05:30:00 IST"
[18] "1980-01-21 05:30:00 IST"
[19] "1980-03-18 05:30:00 IST"
[20] "1980-06-20 05:30:00 IST"
[21] "1980-09-08 05:30:00 IST"
[22] "1981-01-27 05:30:00 IST"
[23] "1981-04-29 05:30:00 IST"
[24] "1981-06-26 05:30:00 IST"
[25] "1981-08-25 05:30:00 IST"
[26] "1982-04-27 05:30:00 IST"
[27] "1982-08-17 05:30:00 IST"
[28] "1982-09-27 05:30:00 IST"
[29] "1983-01-24 05:30:00 IST"
[30] "1983-07-12 05:30:00 IST"
[31] "1984-02-09 05:30:00 IST"
[32] "1984-03-12 05:30:00 IST"
[33] "1984-07-05 05:30:00 IST"
[34] "1984-10-19 05:30:00 IST"
[35] "1985-03-25 05:30:00 IST"
[36] "1985-07-10 05:30:00 IST"
[37] "1985-09-13 05:30:00 IST"
[38] "1985-11-05 05:30:00 IST"
[39] "1986-02-27 05:30:00 IST"
[40] "1986-04-28 05:30:00 IST"
[41] "1986-06-03 05:30:00 IST"
[42] "1986-07-04 05:30:00 IST"
[43] "1986-09-05 05:30:00 IST"
[44] "1987-02-04 05:30:00 IST"
[45] "1987-04-27 05:30:00 IST"
[46] "1987-06-02 05:30:00 IST"
[47] "1987-11-09 05:30:00 IST"
[48] "1988-03-25 05:30:00 IST"
[49] "1988-05-11 05:30:00 IST"
[50] "1989-10-16 05:30:00 IST"
[51] "1990-01-04 05:30:00 IST"
[52] "1990-02-26 05:30:00 IST"
[53] "1990-04-23 05:30:00 IST"
[54] "1990-09-25 05:30:00 IST"
[55] "1991-01-14 05:30:00 IST"
[56] "1991-02-26 05:30:00 IST"
[57] "1991-05-17 05:30:00 IST"
[58] "1991-06-28 05:30:00 IST"
[59] "1991-08-19 05:30:00 IST"
[60] "1992-05-13 05:30:00 IST"
[61] "1992-09-24 05:30:00 IST"
[62] "1993-01-25 05:30:00 IST"
[63] "1993-05-14 05:30:00 IST"
[64] "1993-07-27 05:30:00 IST"
[65] "1993-11-22 05:30:00 IST"
[66] "1994-06-20 05:30:00 IST"
[67] "1994-11-23 05:30:00 IST"
[68] "1995-03-09 05:30:00 IST"
[69] "1995-09-22 05:30:00 IST"
[70] "1995-10-23 05:30:00 IST"
> out <- pot(danish,10)
> plot(out)
https://www.mathclasstutor.online

https://www.mathclasstutor.online

https://www.mathclasstutor.online

https://www.mathclasstutor.online

https://www.mathclasstutor.online

https://www.mathclasstutor.online

https://www.mathclasstutor.online

Make a plot selection (or 0 to exit):

1: plot: Point Process of Exceedances
2: plot: Scatterplot of Gaps
3: plot: Qplot of Gaps
4: plot: ACF of Gaps
5: plot: Scatterplot of Residuals
6: plot: Qplot of Residuals
7: plot: ACF of Residuals
8: plot: Go to GPD Plots

Selection: 1

Make a plot selection (or 0 to exit):

1: plot: Point Process of Exceedances
2: plot: Scatterplot of Gaps
3: plot: Qplot of Gaps
4: plot: ACF of Gaps
5: plot: Scatterplot of Residuals
6: plot: Qplot of Residuals
7: plot: ACF of Residuals
8: plot: Go to GPD Plots

Selection: 2

Make a plot selection (or 0 to exit):

1: plot: Point Process of Exceedances
2: plot: Scatterplot of Gaps
3: plot: Qplot of Gaps
4: plot: ACF of Gaps
5: plot: Scatterplot of Residuals
6: plot: Qplot of Residuals
7: plot: ACF of Residuals
8: plot: Go to GPD Plots

Selection: 3

Make a plot selection (or 0 to exit):

1: plot: Point Process of Exceedances
2: plot: Scatterplot of Gaps
3: plot: Qplot of Gaps
4: plot: ACF of Gaps
5: plot: Scatterplot of Residuals
6: plot: Qplot of Residuals
7: plot: ACF of Residuals
8: plot: Go to GPD Plots

Selection: 4

Make a plot selection (or 0 to exit):

1: plot: Point Process of Exceedances
2: plot: Scatterplot of Gaps
3: plot: Qplot of Gaps
4: plot: ACF of Gaps
5: plot: Scatterplot of Residuals
6: plot: Qplot of Residuals
7: plot: ACF of Residuals
8: plot: Go to GPD Plots

Selection: 5

Make a plot selection (or 0 to exit):

1: plot: Point Process of Exceedances
2: plot: Scatterplot of Gaps
3: plot: Qplot of Gaps
4: plot: ACF of Gaps
5: plot: Scatterplot of Residuals
6: plot: Qplot of Residuals
7: plot: ACF of Residuals
8: plot: Go to GPD Plots

Selection: 6

Make a plot selection (or 0 to exit):

1: plot: Point Process of Exceedances
2: plot: Scatterplot of Gaps
3: plot: Qplot of Gaps
4: plot: ACF of Gaps
5: plot: Scatterplot of Residuals
6: plot: Qplot of Residuals
7: plot: ACF of Residuals
8: plot: Go to GPD Plots

Selection: 7

Make a plot selection (or 0 to exit):

1: plot: Point Process of Exceedances
2: plot: Scatterplot of Gaps
3: plot: Qplot of Gaps
4: plot: ACF of Gaps
5: plot: Scatterplot of Residuals
6: plot: Qplot of Residuals
7: plot: ACF of Residuals
8: plot: Go to GPD Plots

Selection: 8

Make a plot selection (or 0 to exit):

1: plot: Excess Distribution
2: plot: Tail of Underlying Distribution
3: plot: Scatterplot of Residuals
4: plot: QQplot of Residuals

Selection: 0

Make a plot selection (or 0 to exit):

1: plot: Point Process of Exceedances
2: plot: Scatterplot of Gaps
3: plot: Qplot of Gaps
4: plot: ACF of Gaps
5: plot: Scatterplot of Residuals
6: plot: Qplot of Residuals
7: plot: ACF of Residuals
8: plot: Go to GPD Plots

Selection: 0
> out <- gpd(danish, 10)
> shape(danish)
> hill(danish)
> hill(danish, option = "quantile", end = 500, p = 0.999)
> quant(danish, 0.999)
> out <- gpd(danish, 10)
> tp <- tailplot(out)
> gpd.q(tp, 0.999)
 Lower CI  Estimate  Upper CI
 64.66184  94.28956 188.91752
> gpd.sfall(tp, 0.999)
 Lower CI  Estimate  Upper CI
 96.64625 191.36972 394.87555 

Intersection of subring

How to prove the intersection of two subrings is a subring? Let S1 and S2  be two subrings of a ring R. Then S1∩S2 is not empty since ...