cleaning up armax code

2 files changed, 21 insertions, 21 deletions

diff --git a/R/estUtil.R b/R/estUtil.R
index 95a387e..ed54da4 100644
--- a/R/estUtil.R
+++ b/R/estUtil.R
@@ -112,32 +112,31 @@ getcov <- function(sys){
 armaxGrad <- function(theta,e,dots){
   y <- dots[[1]]; u <- dots[[2]]; order <- dots[[3]];
   na <- order[1];nb <- order[2]; nc <- order[3]; nk <- order[4]
-  nb1 <- nb+nk-1 ; n <- max(na,nb1,nc)
-  
-  N <- dim(y)[1]-2*n
+  nb1 <- nb+nk-1 ; n <- max(na,nb1,nc);N <- dim(y)[1]
   
+  l <- list()
   if(is.null(e)){
-    e <- matrix(c(dots[[4]][,],rep(0,n)))
-    eout <- matrix(c(rep(0,n),e[,]))
-  } else{
-    eout <- matrix(c(rep(0,n),e[,]))
-  }
+    e <- dots[[4]]; l$e <- e
+  } 
+  
+  yout <- apply(y,2,padZeros,n=n)
+  uout <- apply(u,2,padZeros,n=n)
+  eout <- apply(e,2,padZeros,n=n)
   
   reg <- function(i) {
     if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1
-    matrix(c(-y[i-1:na,],u[v,],eout[i-1:nc,]))
+    matrix(c(-yout[i-1:na,],uout[v,],eout[i-1:nc,]))
   }
   
   X <- t(sapply(n+1:(N+n),reg))
-  Y <- y[n+1:(N+n),,drop=F]
-  l <- list(X=X,Y=Y,e=e)
+  Y <- yout[n+1:(N+n),,drop=F]
+  fn <- Y-X%*%theta
   
-  if(!is.null(e)){
-    filt1 <- signal::Arma(b=1,a=c(1,theta[(na+nb+1:nc)]))
-    grad <- apply(X,2,signal::filter,filt=filt1) 
-    l$grad <- grad
-  }
+  # Compute Gradient
+  filt1 <- signal::Arma(b=1,a=c(1,theta[(na+nb+1:nc)]))
+  grad <- apply(X,2,signal::filter,filt=filt1) 
   
+  l$grad <- grad[1:N,,drop=F];l$fn <- fn[1:N,,drop=F]
   return(l)
 }
 
@@ -223,4 +222,5 @@ checkInitSys <- function(init_sys){
   }
 }
 
-leftPadZeros <- function(x,n) c(rep(0,n),x)
\ No newline at end of file
+leftPadZeros <- function(x,n) c(rep(0,n),x)
+padZeros <- function(x,n) c(rep(0,n),x,rep(0,n))
\ No newline at end of file
diff --git a/R/estpoly.R b/R/estpoly.R
index 538e7d4..93e50a3 100644
--- a/R/estpoly.R
+++ b/R/estpoly.R
@@ -255,9 +255,9 @@ armax <- function(x,order=c(0,1,1,0),options=optimOptions()){
   if(nc<1) 
     stop("Error: Not an ARMAX model")
   
-  padZeros <- function(x,n) c(rep(0,n),x,rep(0,n))
-  yout <- apply(y,2,padZeros,n=n)
-  uout <- apply(u,2,padZeros,n=n)
+#   padZeros <- function(x,n) c(rep(0,n),x,rep(0,n))
+#   yout <- apply(y,2,padZeros,n=n)
+#   uout <- apply(u,2,padZeros,n=n)
   
   # Initial Parameter Estimates
   mod_arx <- iv(x,c(na,nb,nk)) # fitting ARX model
@@ -266,7 +266,7 @@ armax <- function(x,order=c(0,1,1,0),options=optimOptions()){
   e_init <- matrix(mod_ma$residuals); e_init[is.na(e_init)] <- 0 
   theta0 <- matrix(c(mod_arx$sys$A[-1],mod_arx$sys$B,mod_ma$coef))
   
-  l <- levbmqdt(yout,uout,order,e_init,obj=armaxGrad,
+  l <- levbmqdt(y,u,order,e_init,obj=armaxGrad,
                 theta0=theta0,N=N,opt=options)
   theta <- l$params
   e <- ts(l$residuals,start = start(y),deltat = deltat(y))