# R/0_sim.R In ExtremalDep: Extremal Dependence Models

#### Defines functions simDataMthMrfthfrsimulationsClover

```###
### Generate the random variable from clover distribution
###

simulationsClover = function(n){

R1 = numeric(n);
Theta = numeric(n);
ind = numeric(n);

for(j in 1:n){
U1=runif(1,0,1);
U2=runif(1,0,1);

R1[j] = ((1-U1)^(-2)-1)^(1/2);
f = function(t){2*t/pi+3*sin(4*t)/(10*pi)-U2}
Theta[j] = uniroot(f,lower=0,upper=pi/2)\$root;
ind[j] = 1;
}

Z = numeric(2*n);
dim(Z) = c(n,2);
Z[,1] = R1*cos(Theta);
Z[,2] = R1*sin(Theta); #    Z is simulated from clover distribution
Z }

###
### Generate the random variable from asymmetric distribution
###

# auxiliary functions

fr = function(r){0.30746/(r^(5/12)*(r+1))}
fth = function(th){0.15739/(th^(3/4)*((th+1)^(7/12)))}

Mr = function(r){if(r <= 1){ret = r^(-5/12)}
if(r > 1){ret = r^(-17/12)}
0.30746*ret}

Mth = function(th){if(th <= 1){ret = th^(-3/4)}
if(th > 1){ret = th^(-4/3)}
0.15739*ret}

fr = Vectorize(fr, vectorize.args = "r")
fth = Vectorize(fth, vectorize.args = "th")

Mr = Vectorize(Mr, vectorize.args = "r")
Mth = Vectorize(Mth, vectorize.args = "th")

simData=function(n){

# Step 1: simulate data from m (its df)

Ur = runif(2*n, 0, 1)
Uth = runif(2*n, 0, 1)

small_r = which(Ur <= 5/12)
big_r = which(Ur > 5/12)

small_th = which(Uth <= 4/7)
big_th = which(Uth > 4/7)

Tr = c(1:(2*n))*0
Tth = c(1:(2*n))*0

Tr[small_r] = ((12/5)*Ur[small_r])^(12/7)
Tr[big_r] = ((12/7)*(1-Ur[big_r]))^(-12/5)

Tth[small_th] = ((7/4)*Uth[small_th])^4
Tth[big_th] = ((7/3)*(1-Uth[big_th]))^(-3)

# Step 2: accept if M(T)xU <= f(T)

Umr = runif(2*n, 0, 1)
Umth = runif(2*n, 0, 1)

accept = which(Mr(Tr)*Umr <= fr(Tr) & Mth(Tth)*Umth <= fth(Tth))

L = length(accept)

if(L>=n){
R = c(1:n)*0
TH = c(1:n)*0

R = Tr[accept[1:n]]
TH = Tth[accept[1:n]]

# i=1
# while(i<=n)
# {prod[i] = R[i]*TH[i]; i+1}

X1 = (R/(TH+1))^(1/3)
X2 = ((R*TH)/(TH+1))^(1/4)
ret = list()
ret\$rez = cbind(X1, X2)
ret\$L = L}

if(L<n){print('try again :-)')}
ret}

##  Simulating the data

# in R Console type:
# X = simData(n)
#		     X gives two outputs:
#		     1. n data points (they are matched
#		        here in the way that the point is rejected
#                   if either of 'r' or 'theta' is a bloop)
#		     2. the number of the accepted data

# X = simData(n)\$rez
#		     this gives only the n data points
```

## Try the ExtremalDep package in your browser

Any scripts or data that you put into this service are public.

ExtremalDep documentation built on Aug. 29, 2019, 9:03 a.m.