Skip to content
Prev 935 / 15274 Next

arma model results when exogenouse variables used in, ARMA(p=5, q=(1-6, 19)) (Joe Byers)

All I have a simply script that replicates by error.  I have include a 
summary.fARMA.new function because the digits variables needed 
modification and printCoefmat did not print correctly.  The code is
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Example to replicate ARMA model with exogenous variables error in the 
fARMA.summary functions
#set up enviroment and libraries
library(date); 
library(car);library(tseries);library(Hmisc);library(fBasics);
library(R2HTML);#library(gdata);
library(partsm);library(tseries);library(pear);
library(fBasics);library(fSeries); library(cba);
summary.fARMA.new<-function (object, doplot = TRUE, ...)
{
    digits <- max(5, getOption("digits") - 4)      #added by joe w. byers
    ans = NULL
    x = object
    object = x at fit
    ans$call = object$call
    ans$tsmodel = object$tstitle
    ans$residuals = as.vector(na.omit(object$residuals))
    if (length(ans$residuals) == 0) {
        ans$var = 0
    }
    if (length(ans$residuals) > 0) {
        ans$var = var(ans$residuals)
    }
    ans$sigma2 = object$sigma2
    # the following lines replace to handle exongeous variables and 
masking of paramters
#    tval <- object$coef/object$se.coef
#    prob <- 2 * (1 - pnorm(abs(tval)))
#    ans$coefmat <- cbind(format(object$coef,digits=digits), 
format(object$se.coef,digits=digits),
#        format(tval,digits=digits), prob=format.pval(prob,digits=digits))
#    rownames(ans$coefmat)<-ans$coefmat$Row.names
#    ans$coefmat<-ans$coefmat[2:length(ans$coefmat)]
#    dimnames(ans$coefmat) <- list(names(object$coef), c(" Estimate",
#        " Std. Error", " t value", "Pr(>|t|)"))
#    row.names(ans$coefmat)<-toupper(row.names(ans$coefmat))
    tval<-subset(object$coef,object$mask)/object$se.coef
    prob <- 2 * (1 - pnorm(abs(tval)))
    ans$coefmat<-merge(format(object$coef,digits=digits),
      
cbind(format(object$se.coef,digits=digits),format(tval,digits=digits),format.pval(prob,digits=digits)),
        by.x="row.names",by.y="row.names",sort=FALSE,all.x=TRUE,all.y=FALSE)
    rownames(ans$coefmat)<-toupper(as.vector(ans$coefmat$Row.names)) 
#make the row indexes the parameter names
    ans$coefmat<-ans$coefmat[2:length(ans$coefmat)] # drop first column 
that merge create of row.names
       names(ans$coefmat) <- c(" Estimate"," Std. Error", " t value", 
"Pr(>|t|)")
    if (object$tsmodel == "ar") {
        ans$aic = (object$n.used * (1 + log(2 * pi)) + object$n.used *
            log(ans$var) + 2 * length(object$coef))
    }
    if (object$tsmodel == "arma") {
        ans$aic = (object$n.used * (1 + log(2 * pi)) + object$n.used *
            log(ans$var) + 2 * length(object$coef))
        ans$css = object$css
    }
    if (object$tsmodel == "arima") {
        ans$aic = object$aic
        ans$loglik = object$loglik
    }
    if (object$tsmodel == "fracdiff") {
        doplot = FALSE
    }
    cat("\nTitle:\n ")
    cat(x at title, "\n")
    cat("\nCall:\n ")
    cat(paste(deparse(object$call), sep = "\n", collapse = "\n"),
        "\n", sep = "")
    cat("\nModel:\n ", object$tstitle, "\n", sep = "")
    cat("\nCoefficient(s):\n")
    digits = max(5, getOption("digits") - 4) # changed 4 to 5 by joe w. 
byers
    print.default(format(object$coef, digits = digits), print.gap = 2,
        quote = FALSE)
    digits = max(5, getOption("digits") - 4) # changed 4 to 5 by joe w. 
byers
    if (length(object$residuals) > 2) {
        cat("\nResiduals:\n")
        rq = structure(quantile(ans$residuals), names = c("Min",
            "1Q", "Median", "3Q", "Max"))
        print(rq, digits = digits)
        cat("\nMoments: \n")
        skewness = sum((ans$residuals - 
mean(ans$residuals))^3/sqrt(var(ans$residuals))^3)/length(ans$residuals)
        kurtosis = sum((ans$residuals - 
mean(ans$residuals))^4/var(ans$residuals)^2)/length(ans$residuals) -
            3
        stats = structure(c(skewness, kurtosis), names = c("Skewness",
            "Kurtosis"))
        print(stats, digits = digits)
    }
    cat("\nCoefficient(s):\n")
    signif.stars = getOption("show.signif.stars")
    print(ans$coefmat, digits = digits, signif.stars = signif.stars,
        ...) #changes printCoefmat to print by joe w. byers
    cat("\n")
    if (x at fit$tsmodel == "ar") {
        cat("sigma^2 estimated as:       ", format(object$var,
            digits = digits), "\n")
        cat("AIC Criterion:              ", format(round(object$aic,
            2)), "\n")
    }
    if (x at fit$tsmodel == "arma") {
        cat("sigma^2 estimated as:       ", format(object$sigma2,
            digits = digits), "\n")
        cat("Conditional Sum-of-Squares: ", format(round(object$css,
            digits = 2)), "\n")
    }
    if (x at fit$tsmodel == "arima") {
        cm = object$call$method
        if (is.null(cm) || cm != "CSS")
            cat("sigma^2 estimated as: ", format(object$sigma2,
                digits = digits), "\nlog likelihood:       ",
                format(round(object$loglik, 2)), "\nAIC 
Criterion:        ",
                format(round(object$aic, 2)), "\n", sep = "")
        else cat("sigma^2 estimated as: ", format(object$sigma2,
            digits = digits), "\npart log likelihood:  ", 
format(round(object$loglik,
            2)), "\n", sep = "")
    }
    if (doplot)
        plot.fARMA(x, ...)
    cat("\nDescription:\n ")
    cat(x at description, "\n\n")
    invisible()
}
# generate white noise vector
set.seed(1)           # so you can reproduce the results
v = rnorm(1000,1,1)    # v contains 1000 iid N(1,1) variates
e=rnorm(1000,0,1)
x = cumsum(v)         # x is a random walk with drift = 1
plot.ts(x)            # have a look if you don't believe me
#generate ARMA(1,(1,3)) process with trend
#ar(1)=.5 MA(1)=.05 MA(3)=-.25 trend=.025 mean=1
trend=.025;ma1=.05; ma3=-.25; ar1=.5;
len<-length(v);
data<-ar1*v[3:(len-1)]+ma1*e[3:(len-1)]+ma3*e[1:(len-3)]
plot.ts(data)

# estimate ARMA model using fARMA()
fixed=c(NA,NA,0,NA,NA);#1 ar terms, 3 ma terms, intercept
ts.results<-armaFit(data~arma(1,3),include.mean=TRUE,optim.control=list(maxit=500),fixed=fixed);

#set the graphics to pause since I am using sciViews
par(ask=TRUE);
#produce results with error using summary.fARMA()
summary(ts.results)

#produce results without error using modified summary.fARMA.new()
summary.fARMA.new(ts.results)
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Spencer Graves wrote:
-------------- next part --------------
A non-text attachment was scrubbed...
Name: joe-byers.vcf
Type: text/x-vcard
Size: 295 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20060713/705147a8/attachment.vcf