minor changes

2 files changed, 60 insertions, 74 deletions

diff --git a/R/estUtil.R b/R/estUtil.R
index 8e0e816..5ebcacf 100644
--- a/R/estUtil.R
+++ b/R/estUtil.R
@@ -1,34 +1,3 @@
-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
-  
-  if(is.null(e)){
-    eout <- matrix(rep(0,N+2*n))
-  } else{
-    eout <- matrix(c(rep(0,n),e[,]))
-  }
-  
-  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,]))
-  }
-  
-  X <- t(sapply(n+1:(N+n),reg))
-  Y <- y[n+1:(N+n),,drop=F]
-  l <- list(X=X,Y=Y)
-  
-  if(!is.null(e)){
-    filt1 <- Arma(b=1,a=c(1,theta[(na+nb+1:nc)]))
-    grad <- apply(X,2,filter,filt=filt1) 
-    l$grad <- grad
-  }
-  
-  return(l)
-}
-
 # Implementation of the Levenberg Marquardt Algorithm
 levbmqdt <- function(...,obj,theta0,N,opt){
   dots <- list(...)
@@ -129,4 +98,62 @@ optimOptions <- function(tol=1e-5,maxIter=20,LMinit=2,LMstep=2){
 #' @export
 getcov <- function(sys){
   sys$stats$vcov
+}
+
+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
+  
+  if(is.null(e)){
+    eout <- matrix(rep(0,N+2*n))
+  } else{
+    eout <- matrix(c(rep(0,n),e[,]))
+  }
+  
+  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,]))
+  }
+  
+  X <- t(sapply(n+1:(N+n),reg))
+  Y <- y[n+1:(N+n),,drop=F]
+  l <- list(X=X,Y=Y)
+  
+  if(!is.null(e)){
+    filt1 <- Arma(b=1,a=c(1,theta[(na+nb+1:nc)]))
+    grad <- apply(X,2,filter,filt=filt1) 
+    l$grad <- grad
+  }
+  
+  return(l)
+}
+
+oeGrad <- function(theta,e,dots){
+  # e - Instrument Variable, not necessarily residuals
+  y <- dots[[1]]; u <- dots[[2]]; order <- dots[[3]];
+  nb <- order[1];nf <- order[2]; nk <- order[3];
+  nb1 <- nb+nk-1 ; n <- max(nb1,nf); df <- N - nb - nf
+  
+  N <- dim(y)[1]-n
+  eout <- matrix(c(rep(0,n),e[,]))
+  
+  reg <- function(i) {
+    if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1
+    matrix(c(uout[v,],-eout[i-1:nf,]))
+  }
+  
+  X <- t(sapply(n+1:(N+n),reg))
+  Y <- y[n+1:(N+n),,drop=F]
+  l <- list(X=X,Y=Y)
+  
+  if(!is.null(e)){
+    filt1 <- Arma(b=1,a=c(1,theta[nb+1:nf]))
+    grad <- apply(X,2,filter,filt=filt1)
+    l$grad <- grad
+  }
+  
+  return(l)
 }
\ No newline at end of file
diff --git a/R/estpoly.R b/R/estpoly.R
index b016d82..b74e9db 100644
--- a/R/estpoly.R
+++ b/R/estpoly.R
@@ -357,51 +357,10 @@ oe <- function(x,order=c(1,1,0)){
   
   leftPadZeros <- function(x,n) c(rep(0,n),x)
 
-  reg <- function(i) {
-    if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1
-    matrix(c(uout[v,],-eout[i-1:nf,]))
-  }
-  
-  # Initialize Algorithm
-  i = 0
+  # Initial Guess
   mod_arx <- arx(x,c(nf,nb,nk)) # fitting ARX model
-  iv <- matrix(predict(mod_arx))
-  e <- resid(mod_arx)
-  theta <- c(coef(mod_arx)$B,coef(mod_arx)$A[-1])
-
+  theta0 <- c(coef(mod_arx)$B,coef(mod_arx)$A[-1])
   uout <- apply(u,2,leftPadZeros,n=n)
   
-  tol <- 10^(-5); sumSqRatio <- 1000; lambda <- 1
-  
-  while (sumSqRatio > tol){
-    sumsq0 <- sum(e^2)
-    # Compute gradient
-    eout <- apply(iv,2,leftPadZeros,n=n)
-    X <- t(sapply(n+1:N,reg))
-    filt1 <- Arma(b=1,a=c(1,theta[nb+1:nf]))
-    grad <- apply(X,2,filter,filt=filt1)
-    
-    # Update Parameters
-    H <- 1/N*(t(grad)%*%grad) + lambda*diag(nb+nf)
-    Hinv <- solve(H)
-    theta <- theta + 1/N*Hinv%*%t(grad)%*%e
-    
-    # Update IVs and residuals
-    iv <- X%*%theta; e <- y-iv
-    sumsq <- sum(e^2)
-    
-    sumSqRatio <- abs(sumsq0-sumsq)/sumsq0
-    #     print(sumsq);print(sumSqRatio)
-    i=i+1
-  }
-  # print(sumSqRatio)
-  sigma2 <- sum(e^2)/df
-  vcov <- sigma2 * Hinv
-  
-  model <- idpoly(B = theta[1:nb],F1 = c(1,theta[nb+1:nf]),
-                  ioDelay = nk,Ts=deltat(x))
   
-  estpoly(coefficients = model,vcov = vcov, sigma = sqrt(sigma2),
-          df = df,fitted.values=y-e, residuals=e[,],call=match.call(),
-          input=u)
 }
\ No newline at end of file