diff options
| -rw-r--r-- | NAMESPACE | 1 | ||||
| -rw-r--r-- | R/estUtil.R | 39 | ||||
| -rw-r--r-- | R/estpoly.R | 41 |
3 files changed, 75 insertions, 6 deletions
@@ -31,6 +31,7 @@ export("inputNames<-") export("outputNames<-") export(armax) export(arx) +export(bj) export(compare) export(dataSlice) export(detrend) diff --git a/R/estUtil.R b/R/estUtil.R index 2aadc0e..dc0a00d 100644 --- a/R/estUtil.R +++ b/R/estUtil.R @@ -159,4 +159,43 @@ oeGrad <- function(theta,e,dots){ } return(l) +} + +bjGrad <- function(theta,e,dots){ + y <- dots[[1]]; uout <- dots[[2]]; order <- dots[[3]]; + 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); + N <- dim(y)[1] + + if(is.null(e)){ + zeta <- dots[[4]] + w <- y-zeta + e <- dots[[5]] + } else{ + filt_ts <- signal::Arma(b=c(1,theta[nb+1:nc]), + a=c(1,theta[nb+nc+1:nd])) + w <- signal::filter(filt_ts,e) + zeta <- y-w + } + zetaout <- matrix(c(rep(0,n),zeta[,])) + wout <- matrix(c(rep(0,n),w[,])) + 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:nc,],wout[i-1:nd,],-zetaout[i-1:nf,])) + } + + X <- t(sapply(n+1:N,reg)) + l <- list(X=X,Y=y) + + if(!is.null(e)){ + filt1 <- signal::Arma(b=c(1,theta[nb+nc+1:nd]), + a=c(1,theta[nb+1:nc])) + grad <- apply(X,2,signal::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 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 @@ -359,4 +354,38 @@ oe <- function(x,order=c(1,1,0),options=optimOptions()){ 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) }
\ No newline at end of file |