Brian Smith
2023-Sep-19 17:39 UTC
[R] Could you manually replicate execution of a R function
Hi,
I have trying to replicate a function from rugarch package manually.
Below is the calculation based on the function,
library(rugarch)
data(dji30ret)
spec = ugarchspec(mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
variance.model = list(model = "gjrGARCH"), distribution.model =
"sstd")
fit = ugarchfit(spec, data = dji30ret[1:1000, 1, drop = FALSE])
spec2 = spec
setfixed(spec2)<-as.list(coef(fit))
filt = ugarchfilter(spec2, dji30ret[1001:2500, 1, drop = FALSE], n.old = 1000)
actual = dji30ret[1001:2500,1]
# location+scale invariance allows to use [mu + sigma*q(p,0,1,skew,shape)]
VaR = fitted(filt) + sigma(filt)*qdist("sstd", p=0.05, mu = 0, sigma =
1,
skew = coef(fit)["skew"], shape=coef(fit)["shape"])
alpha = 0.05
conf.level = 0.95
VaRDurTest(0.05, actual, VaR)$b ## 1.019356
Now I want to replicate the calculation manually as below,
VaR.ind = ifelse(actual < VaR, 1, 0)
N = sum(VaR.ind)
TN = length(VaR.ind)
D = diff(which(VaR.ind == 1))
C = rep(0, length(D))
N = length(D)
pweibull = function (D, a, b, survival = FALSE)
{
cdf = 1 - exp(-(a * D)^b)
if (survival)
cdf = 1 - cdf
return(cdf)
}
dweibull = function (D, a, b, log = FALSE)
{
pdf = (b * log(a) + log(b) + (b - 1) * log(D) - (a * D)^b)
if (!log)
pdf = exp(pdf)
return(pdf)
}
likDurationW = function (pars, D, C, N)
{
b = pars[1]
a = ((N - C[1] - C[N])/(sum(D^b)))^(1/b)
lik = C[1] * log(pweibull(D[1], a, b, survival = TRUE)) +
(1 - C[1]) * dweibull(D[1], a, b, log = TRUE) + sum(dweibull(D[2:(N -
1)], a, b, log = TRUE)) + C[N] * log(pweibull(D[N],
a, b, survival = TRUE)) + (1 - C[N]) * dweibull(D[N],
a, b, log = TRUE)
if (!is.finite(lik) || is.nan(lik))
lik = 10000000000
else lik = -lik
return(lik)
}
optim(par = 2, fn = likDurationW, gr = NULL, D = D, C = C, N = N,
method = "L-BFGS-B", lower = 0.001, upper = 10, control = list(trace
0))$par ## 1.029628
I fail to understand why there has to be a difference. Could you
please help to find the reason?
> sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur ... 10.16
Matrix products: default
BLAS:
/Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK:
/Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
locale:
[1] C/UTF-8/C/C/C/C
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] rugarch_1.4-9
loaded via a namespace (and not attached):
[1] Rcpp_1.0.9 mclust_6.0.0
[3] lattice_0.20-45 mvtnorm_1.1-3
[5] zoo_1.8-12 MASS_7.3-58.1
[7] GeneralizedHyperbolic_0.8-4 truncnorm_1.0-9
[9] grid_4.2.2 pracma_2.4.2
[11] KernSmooth_2.23-20 SkewHyperbolic_0.4-0
[13] Matrix_1.5-1 xts_0.12.1
[15] spd_2.0-1 tools_4.2.2
[17] ks_1.14.1 numDeriv_2016.8-1.1
[19] compiler_4.2.2 DistributionUtils_0.6-1
[21] Rsolnp_1.16
Brian Smith
2023-Sep-20 11:42 UTC
[R] Unable to manually replicate execution of a R function
Hi,
** In may earlier post there were some typo, so reposting the same
after correction**
I am trying to replicate a function from rugarch package manually.
Below is the calculation based on the function,
library(rugarch)
data(dji30ret)
spec = ugarchspec(mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
variance.model = list(model = "gjrGARCH"), distribution.model =
"sstd")
fit = ugarchfit(spec, data = dji30ret[1:1000, 1, drop = FALSE])
spec2 = spec
setfixed(spec2)<-as.list(coef(fit))
filt = ugarchfilter(spec2, dji30ret[1001:2500, 1, drop = FALSE], n.old = 1000)
actual = dji30ret[1001:2500,1]
# location+scale invariance allows to use [mu + sigma*q(p,0,1,skew,shape)]
VaR = fitted(filt) + sigma(filt)*qdist("sstd", p=0.05, mu = 0, sigma =
1,
skew = coef(fit)["skew"], shape=coef(fit)["shape"])
alpha = 0.05
conf.level = 0.95
VaRDurTest(0.05, actual, VaR)$b ## 1.019356
Now I want to replicate the calculation manually as below,
VaR.ind = ifelse(actual < VaR, 1, 0)
N = sum(VaR.ind)
TN = length(VaR.ind)
D = diff(which(VaR.ind == 1))
C = rep(0, length(D))
N = length(D)
pweibull = function (D, a, b, survival = FALSE)
{
cdf = 1 - exp(-(a * D)^b)
if (survival)
cdf = 1 - cdf
return(cdf)
}
dweibull = function (D, a, b, log = FALSE)
{
pdf = (b * log(a) + log(b) + (b - 1) * log(D) - (a * D)^b)
if (!log)
pdf = exp(pdf)
return(pdf)
}
likDurationW = function (pars, D, C, N)
{
b = pars[1]
a = ((N - C[1] - C[N])/(sum(D^b)))^(1/b)
lik = C[1] * log(pweibull(D[1], a, b, survival = TRUE)) +
(1 - C[1]) * dweibull(D[1], a, b, log = TRUE) + sum(dweibull(D[2:(N -
1)], a, b, log = TRUE)) + C[N] * log(pweibull(D[N],
a, b, survival = TRUE)) + (1 - C[N]) * dweibull(D[N],
a, b, log = TRUE)
if (!is.finite(lik) || is.nan(lik))
lik = 10000000000
else lik = -lik
return(lik)
}
optim(par = 2, fn = likDurationW, gr = NULL, D = D, C = C, N = N,
method = "L-BFGS-B", lower = 0.001, upper = 10, control = list(trace
0))$par ## 1.029628
I fail to understand why there has to be a difference. Could you
please help to find the reason?
> sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur ... 10.16
Matrix products: default
BLAS:
/Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK:
/Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
locale:
[1] C/UTF-8/C/C/C/C
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] rugarch_1.4-9
loaded via a namespace (and not attached):
[1] Rcpp_1.0.9 mclust_6.0.0
[3] lattice_0.20-45 mvtnorm_1.1-3
[5] zoo_1.8-12 MASS_7.3-58.1
[7] GeneralizedHyperbolic_0.8-4 truncnorm_1.0-9
[9] grid_4.2.2 pracma_2.4.2
[11] KernSmooth_2.23-20 SkewHyperbolic_0.4-0
[13] Matrix_1.5-1 xts_0.12.1
[15] spd_2.0-1 tools_4.2.2
[17] ks_1.14.1 numDeriv_2016.8-1.1
[19] compiler_4.2.2 DistributionUtils_0.6-1
[21] Rsolnp_1.16
Ivan Krylov
2023-Sep-21 07:34 UTC
[R] Could you manually replicate execution of a R function
On Tue, 19 Sep 2023 23:09:18 +0530 Brian Smith <briansmith199312 at gmail.com> wrote:> C = rep(0, length(D)) > N = length(D)In the VaRDurTest function, there's additional code between these two expressions that deals with censoring if head(VaR.ind, 1) == 0 or tail(VaR.ind, 1) == 0. Both of these are true for your input data. -- Best regards, Ivan