[R] farimaSim

[Diethelm Wuertz]

Hansi Weissensteiner wrote:

>Hello!
>
>I installed the fSeries package to get some farima time-series which i tried
>with farimaSim, but unfortunately i got always an error. I tried it this way:
>
>  
>
>>farimaSim(n = 1000, model = list(ar = 0.5,  d = 0.3, ma = 0.1), method="freq")
>>    
>>
>
>Error in farimaSim(n = 1000, model = list(ar = 0.5, d = 0.3, ma = 0.1),  :
> ... used in an incorrect context
>  
>
In the function definition of farimaSim the dots argument is missing.
Replace or overwrite  farimaSim by the following function.
I have to aplogize for any inconveniances caused by this bug.
The bug will be fixed in the next version of Rmetrics.

Diethelm Wuertz



farimaSim =
function(n = 1000, model = list(ar = c(0.5, -0.5), d = 0.3, ma = 0.1),
method = c("freq", "time"), ...)
############  ^^^^^^^^^^   #### Insert the dots!!!!
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Simulates a FARMA Time Series Process
   
    # Note:
    #   Splus-Like argument list
   
    # Example:
    #   armaSim(model = list(ar = c(0.5, -0.5), d = 0.2, ma = 0.1))
    #   armaSim(model = list(d = 0.2, ma = 0))
    #   armaSim(model = list(d = 0.2))
   
    # FUNCTION:
   
    # Settings:
    innov = NULL
    n.start = 100
    start.innov = NULL
    rand.gen = rnorm
   
    # Simulate:
    if (!is.list(model))
        stop("model must be list")
    if (is.null(innov))
        innov = rand.gen(n, ...)
    n = length(innov)
    if (is.null(start.innov))
        start.innov = rand.gen(n, ...)
    n.start = length(start.innov)

    # AR PART:
    p = length(model$ar)
    if (p == 1 && model$ar == 0)
        p = 0
    if (p) {
        minroots = min(Mod(polyroot(c(1, -model$ar))))
        if (minroots <= 1) stop("ar part of model is not stationary")
    }
   
    # MA PART:
    q = length(model$ma)
    if (q == 1 && model$ma == 0)
        q = 0
    if (n.start < p + q)
        stop("burn-in must be as long as ar + ma")
   
    # DIFFERENCING:
    ## if (model$d < 0) stop("d must be positive ")
    dd = length(model$d)   
    if (dd) {
        # FRACDIFF if "dd" is a non-integer value:
        d = model$d
        if (d != round(d) ) {
            TSMODEL = "FRACDIFF"
        } else {
            TSMODEL = "ARIMA" }
    } else {
        d = 0
        TSMODEL = "ARIMA"
    }
   
    # ARMA:
    if (TSMODEL == "ARIMA") {
        stop("d is a short range model")
    }
       
    if (TSMODEL == "FRACDIFF") {
        if (p == 0) model$ar = 0
        if (q == 0) model$ma = 0
        mu = 0
        # Use Fortran Routine from R's contributed fracdiff package:
        # This is a BUILTIN function ...
        x = .Fortran("fdsim", as.integer(n), as.integer(p), as.integer(q),
            as.double(model$ar), as.double(model$ma), as.double(model$d),
            as.double(mu), as.double(rnorm(n + q)), x = double(n + q),
            as.double(.Machine$double.xmin), 
as.double(.Machine$double.xmax),
            as.double(.Machine$double.neg.eps), 
as.double(.Machine$double.eps),
            PACKAGE = "fSeries")$x[1:n]
    }
              
    # Return Value:
    ans = as.ts(x)
    attr(ans, "model") = model
    ans
}



The result will be:

farimaSim(n = 1000, model = list(ar = 0.5,  d = 0.3, ma = 0.1), 
method="freq")
Time Series:
Start = 1
End = 1000
Frequency = 1
   [1]  3.2301357515  2.7050870252  3.4878581467  4.1825601460  5.1286599175
   [6]  5.5220029973  5.1488261085  5.9423142081  3.5078481961  2.1189096844
  [11]  2.5885700846  2.2224732255  3.0389690791  2.8123050401  2.6096020908
  [16]  2.3095398924  1.4900953088  2.6017795027  3.5148157279  3.4317135045
...
 [991] -0.5745174996  0.9022402358 -0.5675451281  1.4723458150  2.2187064082
 [996]  1.6818662030  0.3170217298  0.9290833661  0.5800528928 -1.6796471062
attr(,"model")
attr(,"model")$ar
[1] 0.5

attr(,"model")$d
[1] 0.3

attr(,"model")$ma
[1] 0.1




>Some ideas?
>
>Regards,
>
>    ___
> _ /_|_|   Hansi Weissensteiner
>/o\__/O\=  csae1552 at uibk.ac.at
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>  
>