In some tests I am getting a difference in the 7th or 8th significant digit
depending
on whether I set DUP=T or DUP=F in the .Fortran call below. While this is not a
huge
error it is bigger than I would expect and may be a symptom of another problem.
Any
comments or suggestions? (R1.2.2 on Solaris.)
Paul Gilbert
_______
genD.ARMA <- function(model, data, d=0.01, eps=1e-4, r=6, warn=F){
n <-length(c(output.data(data))) # this has the same length as the
residual
sampleT <-periods(data)
ps <-dim(model$A)[3]
ms <-dim(model$C)[3]
ns <-ps # this could be 1 except z0 is used for TREND and F, Q and R
for
scratch
loc <- match(model$location,
c("A","B","C","t"))
cloc <- match(model$const.location,
c("A","B","C","t"))
if(is.null(ms))
{ms<-0
C<-array(0,c(1,ns,1)) # can't call compiled with 0 length arrays
u <- matrix(0,sampleT,1)
G <- matrix(0,ns,1)
}
else
{C <-model$C
u <- input.data(data)
G <- matrix(0,ns,ms)
}
x <-model$parms
h0 <- abs(d*x)+eps*(x==0.0)
D <- matrix(1e20, n,(length(x)*(length(x) + 3))/2)
Daprox <- matrix(0,n,r)
Hdiag <- matrix(0,n,length(x))
Haprox <- matrix(0,n,r)
storage.mode(D) <-"double"
storage.mode(Daprox)<-"double"
storage.mode(Haprox)<-"double"
storage.mode(Hdiag )<-"double"
D <-.Fortran("gend",
D=D,
as.integer(is.ARMA(model)),
p=as.integer(length(x)),
x=as.double(x),
as.double(h0),
as.integer(n), #6
as.integer((length(x)*(length(x)+ 3))/2), #cols of D
f0=double(n),
r=as.integer(r),
# work space for GEND
Haprox=Haprox, #10
Hdiag=Hdiag,
Daprox=Daprox,
x=double(length(x)),
delta=double(length(x)),
f1=double(n),
f2=double(n),
# work space for ARMAp / KFp
# cov=matrix(1e20,ps,ps),
# pred is f0,f1,f2 passed above
as.integer(ms), # input dimension m #17
as.integer(ps), # output dimension p
as.integer(sampleT),
as.integer(periods(data)),
u=as.double(u),
y=as.double(output.data(data)),
as.integer(loc), #as.character(model$location), #23
as.integer(model$i),
as.integer(model$j),
as.integer(length(model$const)),
const=as.double(model$const),
as.integer(cloc), #as.character(model$const.location),
as.integer(model$const.i ),
as.integer(model$const.j),
# for ARMA models:
as.integer(model$l), #31
as.integer(model$const.l),
as.integer( dim(model$A)[1]),#1+order of A
as.integer( dim(model$B)[1]),#1+order of B
as.integer( dim(C)[1]),#1+order of C
A=model$A,
B=model$B,
C=C,
# for state space models:
as.integer(ns), # state dim. #39
# state=as.double(matrix(0,sampleT,ns)),
# track=as.double(array(0,c(sampleT,ns,ns))),
z0=as.double(rep(0,ps)), # note: this is TREND for ARMA
P0=as.double(diag(0,ps)),
F=as.double(matrix(0,ns,ns)),
G=as.double(G),
H=as.double(matrix(0,ps,ns)),
K=as.double(matrix(0,ns,ps)),
Q=as.double(matrix(0,ns,ns)),
R=as.double(matrix(0,ps,ps)),
gain=as.integer(F), #48
DUP=T)[c("D","p","f0", "x",
"r")]
D$d <- d
D$eps <- eps
D$f0 <- l(model, data, result="pred") - c(output.data(data))
invisible(classed(D,"Darray")) #constructor
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
