From 483bd07eb2b1f83e1e1428007a8678793e4d7e4b Mon Sep 17 00:00:00 2001
From: Suraj Yerramilli
Date: Wed, 10 Feb 2016 13:17:38 +0530
Subject: added the bj routine

---
 R/estpoly.R | 41 +++++++++++++++++++++++++++++++++++------
 1 file changed, 35 insertions(+), 6 deletions(-)

(limited to 'R/estpoly.R')

diff --git a/R/estpoly.R b/R/estpoly.R
index 6d3a13f..a1faff9 100644
--- a/R/estpoly.R
+++ b/R/estpoly.R
@@ -332,7 +332,6 @@ armax <- function(x,order=c(0,1,1,0),options=optimOptions()){
 #' 
 #' @export
 oe <- function(x,order=c(1,1,0),options=optimOptions()){
-  require(signal)
   y <- outputData(x); u <- inputData(x); N <- dim(y)[1]
   nb <- order[1];nf <- order[2]; nk <- order[3];
   nb1 <- nb+nk-1 ; n <- max(nb1,nf); df <- N - nb - nf
@@ -341,13 +340,9 @@ oe <- function(x,order=c(1,1,0),options=optimOptions()){
     stop("Not an OE model")
   
   leftPadZeros <- function(x,n) c(rep(0,n),x)
-
-  # Initial Guess
-  mod_arx <- arx(x,c(nf,nb,nk)) # fitting ARX model
-  iv <- matrix(predict(mod_arx))
-  theta0 <- matrix(c(mod_arx$sys$B,mod_arx$sys$A[-1]))
   uout <- apply(u,2,leftPadZeros,n=n)
   
+  
   l <- levbmqdt(y,uout,order,iv,obj=oeGrad,theta0=theta0,N=N,
                 opt=options)
   theta <- l$params
@@ -356,6 +351,40 @@ oe <- function(x,order=c(1,1,0),options=optimOptions()){
   model <- idpoly(B = theta[1:nb],F1 = c(1,theta[nb+1:nf]),
                   ioDelay = nk,Ts=deltat(x))
   
+  estpoly(sys = model,stats=list(vcov = l$vcov, sigma = l$sigma),
+          fitted.values=y-e,residuals=e,call=match.call(),input=u,
+          options = options,termination = l$termination)
+}
+
+#' @export
+bj <- function(x,order=c(1,1,1,1,0),init_sys=NULL,
+               options=optimOptions()){
+  y <- outputData(x); u <- inputData(x); N <- dim(y)[1]
+  nb <- order[1];nc <- order[2]; nd <- order[3];
+  nf <- order[4]; nk <- order[5];
+  nb1 <- nb+nk-1 ; n <- max(nb1,nc,nd,nf); df <- N-nb-nc-nd-nf
+  
+  # Initial Guess
+  mod_oe <- oe(x,c(nb,nf,nk))
+  v <- resid(mod_oe); zeta <- predict(mod_oe)
+  mod_arma <- arima(v,order=c(nd,0,nc),include.mean = F)
+  theta0 <- matrix(c(mod_oe$sys$B,coef(mod_arma)[nd+1:nc],
+                  coef(mod_arma)[1:nd],mod_oe$sys$F1[-1]))
+  eps <- resid(mod_arma)
+  
+  leftPadZeros <- function(x,n) c(rep(0,n),x)
+  uout <- apply(u,2,leftPadZeros,n=n) 
+  
+  l <- levbmqdt(y,uout,order,zeta,eps,obj=bjGrad,theta0=theta0,N=N,
+                opt=options)
+  theta <- l$params
+  e <- ts(l$residuals,start = start(y),deltat = deltat(y))
+  
+  model <- idpoly(B = theta[1:nb],C=c(1,theta[nb+1:nc]),
+                  D=c(1,theta[nb+nc+1:nd]),
+                  F1 = c(1,theta[nb+nc+nd+1:nf]),
+                  ioDelay = nk,Ts=deltat(x))
+  
   estpoly(sys = model,stats=list(vcov = l$vcov, sigma = l$sigma),
           fitted.values=y-e,residuals=e,call=match.call(),input=u,
           options = options,termination = l$termination)
-- 
cgit