diff options
| -rw-r--r-- | R/estUtil.R | 89 | ||||
| -rw-r--r-- | R/estpoly.R | 45 |
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 |