Jaimie,
It's really quite easy to "open up" the source and look at the code if
you want to make these checks...
Because of the non-differentiability of the absolute value function, an
approximation is used as defined on page 85 of the paper by Hentschel
(the rugarch package follows the exact formulation of the author).
There are 2 errors in your code.
1. You do not use the approximation (error very small)
2. Because delta=lambda, the code uses a special switch:
kdelta=delta*indicator+lambda. In any case, replace delta (which is zero
in the code) with lambda:
Here is the correct code snippet to use (using your notation
and some simplifications for easier reading):
#############################################################
for(i in 2:101)
{
h[i] =
omega+alpha*(h[i-1]^lambda)*(abs(z[i-1]-eta2)-eta1*(z[i-1]-eta2))^lambda+sum(beta*h[i-1]^lambda)
h[i] = h[i]^(1/lambda)
epsilon[i] = h[i]*z[i]
# since its arma(0,0):
y[i] = mu+epsilon[i]
}
h=h[-1]
epsilon=epsilon[-1]
y=y[-1]
all.equal(as.numeric(h),as.numeric(sigma(sim.ruGarch)))
[1] "Mean relative difference: 2.601269e-06"
# That's the approximation error.
#############################################################
# Use this to get exact match:
for(i in 2:101)
{
h[i]= omega + alpha * ( h[i-1]^lambda)*( (sqrt(0.001^2 + ( z[i-1] -
eta2)^2)- eta1* (z[i-1] - eta2))^lambda) + sum(beta*h[i-1]^lambda )
h[i] = h[i]^(1/lambda)
epsilon[i] = h[i]*z[i]
y[i] = mu+epsilon[i]
}
h=h[-1]
epsilon=epsilon[-1]
y=y[-1]
all.equal(as.numeric(h),as.numeric(sigma(sim.ruGarch)))
>TRUE
-Alexios
On 25/02/2014 05:13, Jaimie Villanueva wrote:
> Hi R users
>
> I'm trying to replicate what the fGarch model does when any of the its
> submodels incorporated are specified. I tried to look for help on the web
> but not succeeded.
>
> Up to now, i am able to get almost all but i got problems to replicate the
> result of the aparch submodel.
>
> I'm following the specification of the model set out in
> "Introduction_to_the_rugarch_package.pdf" (equation 26, to be more
> precise) . It seems like i am not correctly setting the constraints of the
> parameters needed. Based on the document, to get the aparch submodel from
> the general equation of the fGarch model, the following constraints have to
> be set:
>
> - delta = lambda
> - eta2 = 0
> - abs(eta1) <= 1
>
> It seems like i'm missing someone of them.
>
> This is what i did.
>
> library(rugarch)
>
> data(sp500ret)
>
> seed=1
>
> ##################################################
> # SIMULATION WITH THE RUGARCH BUILT-IN FUNCTIONS
>
> spec=ugarchspec(variance.model = list(model = "fGARCH", submodel =
> "APARCH",
> garchOrder = c(1,1)),
> mean.model = list(armaOrder = c(0,0), include.mean = TRUE),
> distribution.model = "norm")
>
> model=ugarchfit(spec,as.xts(sp500ret),solver = "hybrid")
>
> sim.ruGarch = ugarchsim(model, n.sim = 100,n.start = 0, m.sim = 1,
> startMethod = 'sample',rseed = seed)
>
>
> ##################################################
> ### MY ATTEMPT TO MANUALLY REPLICATE IT.
>
>
> # STORE PARAMETERS TO USE THEM LATER ON
>
> ParIndx=model@model$pidx
> Param=model@model$pars
>
> mu=Param[ParIndx["mu",1]:ParIndx["mu",2],1]
> ar=Param[ParIndx["ar",1]:ParIndx["ar",2],1]
> orden.ar=length(ar)
> ma=Param[ParIndx["ma",1]:ParIndx["ma",2],1]
> orden.ma=length(ma)
> arfima=Param[ParIndx["arfima",1]:ParIndx["arfima",2],1]
> archm=Param[ParIndx["archm",1]:ParIndx["archm",2],1]
> mxreg=Param[ParIndx["mxreg",1]:ParIndx["mxreg",2],1]
> omega=Param[ParIndx["omega",1]:ParIndx["omega",2],1]
> alpha=Param[ParIndx["alpha",1]:ParIndx["alpha",2],1]
> orden.alpha=length(alpha)
> beta=Param[ParIndx["beta",1]:ParIndx["beta",2],1]
> orden.beta=length(beta)
> gamma=Param[ParIndx["gamma",1]:ParIndx["gamma",2],1]
> eta1=Param[ParIndx["eta1",1]:ParIndx["eta1",2],1]
> eta2=Param[ParIndx["eta2",1]:ParIndx["eta2",2],1]
> delta=Param[ParIndx["delta",1]:ParIndx["delta",2],1]
> lambda=Param[ParIndx["lambda",1]:ParIndx["lambda",2],1]
> vxreg=Param[ParIndx["vxreg",1]:ParIndx["vxreg",2],1]
> skew=Param[ParIndx["skew",1]:ParIndx["skew",2],1]
> shape=Param[ParIndx["shape",1]:ParIndx["shape",2],1]
> ghlambda=Param[ParIndx["ghlambda",1]:ParIndx["ghlambda",2],1]
>
> # SET INITIAL VALUES.
>
> h=array(0,c(101))
> epsilon=array(0,c(101))
> z=array(0,c(101))
> y=array(0,c(101))
>
> y[1]=tail(model@model$modeldata$data,1)
> h[1]=tail(model@fit$sigma,1)
> epsilon[1]=tail(model@fit$residuals,1)
> z[1]=epsilon[1]/h[1]
>
>
> set.seed(seed)
>
> z[2:101]=rdist(distribution = "norm", n=100, mu = 0, sigma = 1,
> lambda = ghlambda, skew = skew,shape = shape)
>
> for(i in 2:101)
> {
>
> h[i] = omega+
>
>
> sum(alpha*h[i-(1:orden.alpha)]^lambda*((abs(z[i-(1:orden.alpha)]-eta2)-eta1*(z[i-(1:orden.alpha)]-eta2))^delta))+
> sum(beta*h[i-(1:orden.beta)]^lambda)
>
> h[i] = h[i]^(1/lambda)
>
> epsilon[i] = h[i]*z[i]
>
> y[i] = mu+sum(ar*(y[i-(1:orden.ar)]-mu))+sum(ma*epsilon[i-(1:orden.ma
> )])+epsilon[i]
>
> }
>
> h=h[-1]
> epsilon=epsilon[-1]
> y=y[-1]
>
>
> After running this, this is what i get:
>
>> head(data.frame(round(h,7),round(sigma(sim.ruGarch),7)))
> round.h..7. round.sigma.sim.ruGarch...7.
> T+1 0.0247601 0.0260287
> T+2 0.0247232 0.0262280
> T+3 0.0246867 0.0246775
> T+4 0.0246505 0.0255587
> T+5 0.0246145 0.0243801
> T+6 0.0245789 0.0229765
>
>
> Could anyone help me to figure out what i am doing wrong?
>
> Any help would be strongly appreciated.
>
> Cheers!!
>
_______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should
go.
