74 files changed, 4467 insertions, 2 deletions

diff --git a/.Rbuildignore b/.Rbuildignore
new file mode 100644
index 0000000..91114bf
--- /dev/null
+++ b/.Rbuildignore
@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..807ea25
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+.Rproj.user
+.Rhistory
+.RData
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..8ec883a
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,13 @@
+Package: sysid
+Type: Package
+Title: System Identification in R
+Version: 1.0.4
+Date: 2017-01-06
+Author: Suraj Yerramilli, Arun Tangirala
+Maintainer: Suraj Yerramilli <surajyerramilli@gmail.com>
+Description: Provides functions for constructing mathematical models of dynamical systems from measured input-output data. 
+License: GPL-3
+Depends: R (>= 3.1)
+Imports:
+    signal,tframe, ggplot2 (>= 2.1.0), reshape2, polynom, bitops, zoo
+RoxygenNote: 5.0.1
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..2e5a9bf
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,96 @@
+# Generated by roxygen2: do not edit by hand
+
+S3method("%=%",lbunch)
+S3method("inputNames<-",idframe)
+S3method("outputNames<-",idframe)
+S3method(inputData,default)
+S3method(inputData,idframe)
+S3method(inputNames,default)
+S3method(inputNames,idframe)
+S3method(outputData,default)
+S3method(outputData,idframe)
+S3method(outputNames,default)
+S3method(outputNames,idframe)
+S3method(plot,estpoly)
+S3method(plot,idframe)
+S3method(plot,idfrd)
+S3method(predict,estpoly)
+S3method(print,estpoly)
+S3method(print,idpoly)
+S3method(print,summary.estpoly)
+S3method(print,summary.idframe)
+S3method(sim,default)
+S3method(sim,idpoly)
+S3method(summary,estpoly)
+S3method(summary,idframe)
+export("%=%")
+export("inputNames<-")
+export("outputNames<-")
+export(armax)
+export(arx)
+export(bj)
+export(compare)
+export(dataSlice)
+export(deltat)
+export(detrend)
+export(estpoly)
+export(etfe)
+export(fitch)
+export(frequency)
+export(g)
+export(getcov)
+export(idframe)
+export(idfrd)
+export(idinput)
+export(idpoly)
+export(impulseest)
+export(impulseplot)
+export(inputData)
+export(inputNames)
+export(iv)
+export(iv4)
+export(misdata)
+export(nInputSeries)
+export(nOutputSeries)
+export(oe)
+export(optimOptions)
+export(outputData)
+export(outputNames)
+export(rarx)
+export(read.idframe)
+export(read.table.idframe)
+export(residplot)
+export(sim)
+export(spa)
+export(step)
+export(time)
+export(trInfo)
+import(bitops)
+import(ggplot2)
+import(polynom)
+import(reshape2)
+import(signal)
+import(tframe)
+importFrom(graphics,par)
+importFrom(graphics,plot)
+importFrom(grid,grid.layout)
+importFrom(grid,grid.newpage)
+importFrom(grid,pushViewport)
+importFrom(grid,viewport)
+importFrom(stats,acf)
+importFrom(stats,ar)
+importFrom(stats,arima)
+importFrom(stats,ccf)
+importFrom(stats,coef)
+importFrom(stats,fitted)
+importFrom(stats,lm)
+importFrom(stats,mvfft)
+importFrom(stats,predict)
+importFrom(stats,resid)
+importFrom(stats,rnorm)
+importFrom(stats,start)
+importFrom(stats,ts)
+importFrom(stats,window)
+importFrom(utils,read.table)
+importFrom(utils,tail)
+importFrom(zoo,na.approx)
diff --git a/R/data.R b/R/data.R
new file mode 100644
index 0000000..89b65ea
--- /dev/null
+++ b/R/data.R
@@ -0,0 +1,117 @@
+#' Data simulated from an ARMAX model 
+#' 
+#' This dataset contains 2555 samples simulated from the following ARMAX model:
+#' \deqn{
+#'  y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + 
+#'  \frac{1-0.3q^{-1}}{1 - 0.5q^{-1}} e[k]
+#' }
+#' 
+#' The model is simulated with a 2555 samples long full-band PRBS input.
+#' The noise variance is set to 0.1
+#' 
+#' @format an \code{idframe} object with 2555 samples, one input and one
+#' output
+#' 
+"armaxsim"             
+
+#' Data simulated from an ARX model 
+#' 
+#' This dataset contains 2555 samples simulated from the following ARX model:
+#' \deqn{
+#'  y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + 
+#'  \frac{1}{1 - 0.5q^{-1}} e[k]
+#' }
+#' 
+#' The model is simulated with a 2555 samples long full-band PRBS input.
+#' The noise variance is set to 0.1
+#' 
+#' @format an \code{idframe} object with 2555 samples, one input and one
+#' output
+#'
+"arxsim"               
+
+#' Data simulated from an BJ model 
+#' 
+#' This dataset contains 2046 samples simulated from the following BJ model:
+#' \deqn{
+#'  y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + 
+#'  \frac{1+0.2q^{-1}}{1 - 0.3q^{-1}} e[k]
+#' }
+#' 
+#' The model is simulated with a 2046 samples long full-band PRBS input.
+#' The noise variance is set to 0.1
+#' 
+#' @format an \code{idframe} object with 2046 samples, one input and one
+#' output
+#'
+"bjsim"                
+
+#' Continuous stirred tank reactor data (idframe)
+#' 
+#' The Process is a model of a Continuous Stirring Tank Reactor, 
+#' where the reaction is exothermic and the concentration is 
+#' controlled by regulating the coolant flow.
+#' \cr
+#' 
+#' Inputs: q, Coolant Flow l/min
+#' Outputs:
+#' \describe{
+#' \item{Ca}{Concentration mol/l}
+#' \item{T}{Temperature Kelvin}}
+#' 
+#' @format an \code{idframe} object with 7500 samples, one input and two
+#' outputs
+"cstr"                 
+
+#' Continuous stirred tank reactor data (data.frame)
+#' 
+#' The Process is a model of a Continuous Stirring Tank Reactor, 
+#' where the reaction is exothermic and the concentration is 
+#' controlled by regulating the coolant flow.
+#' \cr
+#' 
+#' Inputs: q, Coolant Flow l/min
+#' Outputs:
+#' \describe{
+#' \item{Ca}{Concentration mol/l}
+#' \item{T}{Temperature Kelvin}}
+#' 
+#' @format an \code{data.frame} object with 7500 rows and three columns:
+#' q, Ca and T
+#' 
+#' @source \url{ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/cstr.dat.gz}
+"cstrData"             
+
+#' Continuous stirred tank reactor data with missing values
+#' 
+#' This dataset is derived from the \code{cstr} dataset with few samples
+#' containing missing values, in one or all variables. It is used to 
+#' demonstrate the capabilities of the \code{misdata} routine.
+#' 
+#' @format an \code{idframe} object with 7500 samples, one input and two
+#' outputs
+#' 
+#' @seealso \code{\link{cstr}}, \code{\link{misdata}}
+"cstr_mis"             
+
+#' Frequency response data
+#' 
+#' This dataset contains frequency response data of an unknown SISO system.
+#' 
+#' @format an \code{idfrd} object with response at 128 frequency points 
+"frd"                  
+
+#' Data simulated from an OE model 
+#' 
+#' This dataset contains 2555 samples simulated from the following OE model:
+#' \deqn{
+#'  y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + e[k]
+#' }
+#' 
+#' The model is simulated with a 2555 samples long full-band PRBS input.
+#' The noise variance is set to 0.1
+#' 
+#' @format an \code{idframe} object with 2555 samples, one input and one
+#' output
+#'
+"oesim"
\ No newline at end of file
diff --git a/R/estUtil.R b/R/estUtil.R
new file mode 100644
index 0000000..084760c
--- /dev/null
+++ b/R/estUtil.R
@@ -0,0 +1,232 @@
+# Implementation of the Levenberg Marquardt Algorithm
+levbmqdt <- function(...,obj,theta0,N,opt){
+  dots <- list(...)
+  
+  # Optimization Parameters
+  tol <- opt$tol; maxIter <- opt$maxIter
+  d <- opt$adv$LMinit; mu <- opt$adv$LMstep
+  
+  df <- N - dim(theta0)[1]
+  
+  # Initialize Algorithm
+  i <- 0
+  l <- obj(theta=theta0,e=NULL,dots)
+  e <- l$e; grad <- l$grad
+  sumsq0 <- sum(e^2)/df
+  
+  # variable to count the number of times objective function is called
+  countObj <- 0
+  sumSqDiff <- 9E-3*sumsq0
+  
+  repeat{
+    i=i+1
+    
+    g <- 1/df*t(grad)%*%e
+    termPar <- norm(g,"2")
+    
+    repeat{
+      # Update Parameters
+      H <- 1/df*t(grad)%*%grad + d*diag(dim(theta0)[1])
+      Hinv <- solve(H);
+      
+      theta <- theta0 + Hinv%*%g
+      
+      # Evaulate sum square error
+      l <- obj(theta,e,dots)
+      sumsq <- sum(l$fn^2)/df
+      sumSqDiff <- sumsq0-sumsq
+      countObj <- countObj + 1
+      
+      if(termPar < tol) break
+      # no major improvement
+      if(abs(sumSqDiff) < 0.01*sumsq0) break
+      
+      # If sum square error with the updated parameters is less than the 
+      # previous one, the updated parameters become the current parameters
+      # and the damping coefficient is reduced by a factor of mu
+      if(sumSqDiff > 0){
+        d <- d/mu
+        theta0 <- theta
+        sumsq0 <- sumsq
+        e <- l$fn; grad <- l$grad
+        break
+      } else{ # increase damping coefficient by a factor of mu
+        d <- d*mu
+      } 
+    }
+    if(termPar < tol) {
+      WhyStop <- "Tolerance"
+      break
+    }
+    if(abs(sumSqDiff) < 0.01*sumsq0){
+      WhyStop <- "No significant change"
+      break
+    } 
+    if(i == maxIter){
+      WhyStop <- "Maximum Iteration Limit"
+      break
+    } 
+  }
+  # theta <- theta0
+  sigma2 <- sum(e^2)/df
+  vcov <- 1/df*Hinv*sigma2
+  
+  list(params=theta,residuals=e,vcov=vcov,sigma = sqrt(sigma2),
+       termination=list(WhyStop=WhyStop,iter=i,FcnCount = countObj,
+                        CostFcn=sumsq0))
+}
+
+#' Create optimization options
+#' 
+#' Specify optimization options that are to be passed to the 
+#' numerical estimation routines
+#' 
+#' @param tol Minimum 2-norm of the gradient (Default: \code{1e-2})
+#' @param maxIter Maximum number of iterations to be performed
+#' @param LMinit Starting value of search-direction length 
+#' in the Levenberg-Marquardt method (Default: \code{0.01})
+#' @param LMstep Size of the Levenberg-Marquardt step (Default: \code{2})
+#' @param display Argument whether to display iteration details or not
+#' (Default: \code{"off"})
+#' 
+#' @export
+optimOptions <- function(tol=1e-2,maxIter=20,LMinit=0.01,LMstep=2,
+                         display=c("off","on")[1]){
+  return(list(tol=tol,maxIter= maxIter, 
+              adv= list(LMinit=LMinit,LMstep=LMstep),display=display))
+}
+
+#' Parameter covariance of the identified model
+#' 
+#' Obtain the parameter covariance matrix of the linear, identified 
+#' parametric model
+#' 
+#' @param sys a linear, identified parametric model
+#' 
+#' @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]
+  
+  l <- list()
+  if(is.null(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(-yout[i-1:na,],uout[v,],eout[i-1:nc,]))
+  }
+  
+  X <- t(sapply(n+1:(N+n),reg))
+  Y <- yout[n+1:(N+n),,drop=F]
+  fn <- Y-X%*%theta
+  
+  # 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)
+}
+
+oeGrad <- function(theta,e,dots){
+  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)
+  N <- dim(y)[1]
+  l <- list()
+  if(is.null(e)){
+    iv <- dots[[4]]
+    fn <- y-iv; l$e <- fn
+  } else{
+    iv <- y-e
+  }
+  uout <- apply(u,2,leftPadZeros,n=n)
+  ivout <- apply(iv,2,leftPadZeros,n=n)
+  
+  reg <- function(i) {
+    if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1
+    matrix(c(uout[v,],-ivout[i-1:nf,]))
+  }
+  # Compute new regressor matrix and residuals
+  X <- t(sapply(n+1:N,reg))
+  fn <- y-X%*%theta
+  
+  # Compute gradient
+  filt1 <- signal::Arma(b=1,a=c(1,theta[nb+1:nf,]))
+  grad <- apply(X,2,signal::filter,filt=filt1)
+  
+  l$fn <- fn; l$grad<-grad
+  return(l)
+}
+
+bjGrad <- function(theta,e,dots){
+  y <- dots[[1]]; u <- 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]
+
+  l <- list()    
+  if(is.null(e)){
+    zeta <- dots[[4]]
+    w <- y-zeta
+    e <- dots[[5]]; l$e <- e
+  } else{
+    filt_ts <- signal::Arma(b=c(1,theta[nb+1:nc]),
+                            a=c(1,theta[nb+nc+1:nd]))
+    w <- matrix(signal::filter(filt_ts,e))
+    zeta <- y-w
+  }
+  
+  uout <- apply(u,2,leftPadZeros,n=n)
+  zetaout <- apply(zeta,2,leftPadZeros,n=n)
+  eout <- apply(e,2,leftPadZeros,n=n)
+  wout <- apply(w,2,leftPadZeros,n=n)
+  
+  reg <- function(i) {
+    if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1
+    ereg <- if(nc==0) NULL else eout[i-1:nc,]
+    matrix(c(uout[v,],ereg,wout[i-1:nd,],-zetaout[i-1:nf,]))
+  }
+  
+  # Compute new regressor matrix and residuals
+  X <- t(sapply(n+1:N,reg))
+  fn <- y-X%*%theta
+  
+  # Computing gradient
+  C_params <- if(nc==0) NULL else theta[nb+1:nc] 
+  den <- as.numeric(polynom::polynomial(c(1,C_params))*
+                        polynom::polynomial(c(1,theta[nb+nc+nd+1:nf])))
+  filt1 <- signal::Arma(b=c(1,theta[nb+nc+1:nd]),
+                          a=den)
+  grad <- apply(X,2,signal::filter,filt=filt1)
+  
+  l$fn <- fn; l$grad <- grad
+  return(l)
+}
+
+checkInitSys <- function(init_sys){
+  z <- strsplit(toString(sys.call(which=-1)),split = ",")[[1]][1]
+  if(init_sys$type!=z){
+    errMes <- paste("An idpoly model of",toupper(z),"structure expected for the",z,"command.")
+    stop(errMes)
+  }
+}
+
+leftPadZeros <- function(x,n) c(rep(0,n),x)
+padZeros <- function(x,n) c(rep(0,n),x,rep(0,n))
+
+integfilter <- function(x){
+  as.numeric(stats::filter(x,filter=c(1,-1),"convolution",sides = 1,
+                           circular = T))
+}
\ No newline at end of file
diff --git a/R/estpoly.R b/R/estpoly.R
new file mode 100644
index 0000000..e37ad1d
--- /dev/null
+++ b/R/estpoly.R
@@ -0,0 +1,636 @@
+#' Estimated polynomial object
+#' 
+#' Estimated discrete-time polynomial model returned from an estimation
+#' routine.
+#' 
+#' @param sys an \code{idpoly} object containing the estimated polynomial
+#' coefficients
+#' @param fitted.values 1-step ahead predictions on the training dataset
+#' @param residuals 1-step ahead prediction errors
+#' @param options optimization specification ser used (applicable for non-linear least
+#' squares)
+#' @param call the matched call
+#' @param stats a list containing estimation statistics
+#' @param termination termination criteria for optimization
+#' @param input input signal of the training data-set
+#' 
+#' @details 
+#' Do not use \code{estpoly} for directly specifing an input-output polynomial model.
+#' \code{\link{idpoly}} is to be used instead
+#' 
+#' @export
+estpoly <- function(sys,fitted.values,residuals,options=NULL,
+                    call,stats,termination=NULL,input){
+  out <- list(sys=sys,fitted.values=fitted.values,
+              residuals=residuals,input=input,call=call,
+              stats=stats,options=options,termination=termination)
+  class(out) <- "estpoly"
+  out
+}
+
+
+#' @export
+print.estpoly <- function(x,...){
+  print(summary(x),...)
+}
+
+#' @export
+summary.estpoly <- function(object,...)
+{
+  model <- object$sys
+  coefs <- params(model)
+  
+  se <- sqrt(diag(getcov(object)))
+  params <- data.frame(Estimated=coefs,se=se)
+  
+  report <- list(fit=fitch(object),params=params)
+  res <- list(model=model,report=report)
+  class(res) <- "summary.estpoly"
+  res
+}
+
+#' Fit Characteristics
+#' 
+#' Returns quantitative assessment of the estimated model as a list
+#' 
+#' @param x the estimated model
+#' 
+#' @return 
+#' A list containing the following elements
+#' 
+#' \item{MSE}{Mean Square Error measure of how well the response of the model fits
+#' the estimation data}
+#' \item{FPE}{Final Prediction Error}
+#' \item{FitPer}{Normalized root mean squared error (NRMSE) measure of how well the 
+#' response of the model fits the estimation data, expressed as a percentage.}
+#' \item{AIC}{Raw Akaike Information Citeria (AIC) measure of model quality}
+#' \item{AICc}{Small sample-size corrected AIC}
+#' \item{nAIC}{Normalized AIC}
+#' \item{BIC}{Bayesian Information Criteria (BIC)}
+#' 
+#' @export
+fitch <- function(x){
+  y <- fitted(x) + resid(x)
+  ek <- as.matrix(resid(x))
+  N <- nrow(ek); np <- length(params(x$sys))
+  
+  # fit characteristics
+  mse <- det(t(ek)%*%ek)/N
+  fpe <- mse*(1+np/N)/(1-np/N)
+  nrmse <- 1 - sqrt(sum(ek^2))/sqrt(sum((y-mean(y))^2))
+  AIC <- N*log(mse) + 2*np + N*dim(matrix(y))[2]*(log(2*pi)+1)
+  AICc <- AIC*2*np*(np+1)/(N-np-1)
+  nAIC <- log(mse) + 2*np/N
+  BIC <- N*log(mse) + N*dim(matrix(y))[2]*(log(2*pi)+1) + np*log(N)
+  
+  list(MSE=mse,FPE=fpe,FitPer = nrmse*100,AIC=AIC,AICc=AICc,nAIC=nAIC,BIC=BIC)
+}
+
+#' @export
+print.summary.estpoly <- function(x,digits=4,...){
+  print(x$model,se=x$report$params[,2],dig=digits)
+  cat("\n Fit Characteristics \n")
+  print(data.frame(x$report$fit),digits=digits)
+}
+
+#' @import ggplot2
+#' @export
+plot.estpoly <- function(x,newdata=NULL,...){
+  
+  if(is.null(newdata)){
+    ypred <- ts(fitted(x),names="Predicted")
+    yact <- ts(fitted(x) + resid(x),names="Actual")
+    time <- time(x$input)
+    titstr <- "Predictions of Model on Training Set"
+  } else{  
+    if(class(newdata)!="idframe") stop("Only idframe objects allowed")
+    ypred <- predict(x,newdata)
+    yact <- outputData(newdata)[,1]
+    time <- time(newdata)
+    titstr <- "Predictions of Model on Test Set"
+  }
+  df <- data.frame(Predicted=ypred,Actual=yact,Time=time)
+  with(df,ggplot(df, aes(x = Actual,y=Predicted)) +  ggtitle(titstr) +
+    geom_abline(intercept=0,slope=1,colour="#D55E00") +  geom_point())
+}
+
+#' Plot residual characteristics
+#' 
+#' Computes the 1-step ahead prediction errors (residuals) for an estimated polynomial
+#' model, and plots auto-correlation of the residuals and the 
+#' cross-correlation of the residuals with the input signals.
+#' 
+#' @param model estimated polynomial model
+#' @param newdata an optional dataset on which predictions are to be computed. If
+#' not supplied, predictions are computed on the training dataset.
+#' @export
+residplot <- function(model,newdata=NULL){
+  if(is.null(newdata)){
+    e <- resid(model); u <- model$input
+  } else{
+    if(class(newdata)!="idframe") stop("Only idframe objects allowed")
+    e <- newdata$output[,1] - predict(model,newdata)[,1]
+    u <- newdata$input
+  }
+  e <- matrix(e)
+  acorr <- acf(e[,],plot = F); ccorr <- ccf(u[,1],e[,],plot = F)
+  par(mfrow=c(2,1),mar=c(3,4,3,2))
+  plot(acorr,ci=0.99,main="ACF of residuals")
+  plot(ccorr,ci=0.99,main="CCF between the input and residuals",ylab="CCF")
+}
+
+#' Estimate ARX Models
+#' 
+#' Fit an ARX model of the specified order given the input-output data 
+#' 
+#' @param x an object of class \code{idframe}
+#' @param order Specification of the orders: the three integer components 
+#' (na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) and 
+#' the input-output delay
+#' @param lambda Regularization parameter(Default=\code{0.1}) 
+#' @param intNoise Logical variable indicating whether to add integrators in
+#' the noise channel (Default=\code{FALSE})
+#' @param fixed list containing fixed parameters. If supplied, only \code{NA} entries 
+#' will be varied. Specified as a list of two vectors, each containing the parameters
+#' of polynomials A and B respectively.
+#' 
+#' @details
+#' SISO ARX models are of the form 
+#' \deqn{
+#'    y[k] + a_1 y[k-1] + \ldots + a_{na} y[k-na] = b_{nk} u[k-nk] + 
+#'    \ldots + b_{nk+nb} u[k-nk-nb] + e[k] 
+#' }
+#' The function estimates the coefficients using linear least squares (with
+#' regularization).
+#' \cr
+#' The data is expected to have no offsets or trends. They can be removed 
+#' using the \code{\link{detrend}} function. 
+#' \cr
+#' To estimate finite impulse response(\code{FIR}) models, specify the first
+#' order to be zero.
+#' 
+#' @return
+#' An object of class \code{estpoly} containing the following elements:
+#'  \item{sys}{an \code{idpoly} object containing the 
+#'    fitted ARX coefficients}
+#'  \item{fitted.values}{the predicted response}
+#'  \item{residuals}{the residuals}
+#'  \item{input}{the input data used}
+#'  \item{call}{the matched call}
+#'  \item{stats}{A list containing the following fields: \cr
+#'      \code{vcov} - the covariance matrix of the fitted coefficients \cr
+#'      \code{sigma} - the standard deviation of the innovations\cr
+#'      \code{df} - the residual degrees of freedom}
+#' 
+#' 
+#' @references
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+#' Theory and Practice}, CRC Press, Boca Raton. Section 21.6.1
+#' 
+#' Lennart Ljung (1999), \emph{System Identification: Theory for the User}, 
+#' 2nd Edition, Prentice Hall, New York. Section 10.1
+#' 
+#' @examples
+#' data(arxsim)
+#' mod_arx <- arx(arxsim,c(1,2,2))
+#' mod_arx
+#' plot(mod_arx) # plot the predicted and actual responses
+#' 
+#' @export
+arx <- function(x,order=c(1,1,1),lambda=0.1,intNoise=FALSE,
+                fixed=NULL){
+  y <- outputData(x); u <- inputData(x)
+  if(intNoise){
+    y <- apply(y,2,diff)
+    u <- apply(u,2,diff)
+  }
+  N <- dim(y)[1]
+  na <- order[1];nb <- order[2]; nk <- order[3]
+  nb1 <- nb+nk-1 ; n <- max(na,nb1); df <- N-na-nb
+  
+  yout <- apply(y,2,padZeros,n=n);
+  uout <- apply(u,2,padZeros,n=n);
+  
+  fixedflag = is.null(fixed)
+  uindex <- nk:nb1
+  if(na!=0) yindex <- 1:na
+
+  if(!fixedflag){
+    # checking for correct specification of fixed parameters
+    fixedA <- NULL;fixedB <- NULL
+    g(fixedA,fixedB) %=% lapply(fixed,length)
+    if(fixedA != na || fixedB != nb)
+      stop("Number of parameters incorrectly specified in 'fixed'")
+    
+    fixedpars <- unlist(fixed)
+    fixedpars <- fixedpars[!is.na(fixedpars)]
+    df <- df + length(fixedpars)
+    
+    fixedpos_B <- which(!is.na(fixed[[2]]))
+    uindex <- uindex[!uindex %in% (nk+fixedpos_B-1)]
+    fixedpos_A <- numeric(0)
+    if(na!=0){
+      fixedpos_A <- which(!is.na(fixed[[1]]))
+      yindex <- yindex[!yindex %in% fixedpos_A]
+    }
+  }
+  
+  reg <- function(i) {
+    # regressor
+    temp <- numeric(0)
+    if(na!=0) temp <- c(temp,-yout[i-yindex,])
+    phi <- t(c(temp,uout[i-uindex,]))
+    l <- list(phi=phi)
+    # fixed regressors
+    if(!fixedflag){
+      temp <- numeric(0)
+      if(length(fixedpos_A)!=0){
+        temp <- c(temp,-yout[i-fixedpos_A,])
+      } 
+      if(length(fixedpos_B)!=0){
+        temp <- c(temp,uout[i-(nk+fixedpos_B),])
+      }
+      l$fixed <- t(temp)
+    }
+    return(l)
+  }
+  
+  temp <- lapply(n+1:(N+n),reg)
+  X <- do.call(rbind,lapply(temp, function(x) x[[1]]))
+  Y <- yout[n+1:(N+n),,drop=F]
+  fixedY <- matrix(rep(0,nrow(Y)))
+  if(!fixedflag){
+    fixedreg <- do.call(rbind,lapply(temp, function(x) x[[2]]))
+    fixedY <- fixedreg%*%fixedpars
+    Y <- Y - fixedY
+  }
+  
+  # lambda <- 0.1
+  inner <- t(X)%*%X + lambda*diag(dim(X)[2])
+  innerinv <- solve(inner)
+  pinv <- innerinv%*% t(X)
+  coef <- pinv%*%Y
+  
+  eps <- X%*%coef
+  sigma2 <- sum(eps^2)/(df+n)
+  vcov <- sigma2 * innerinv
+  fit <- (X%*%coef+fixedY)[1:N,,drop=F]
+  if(intNoise) fit <- apply(fit,2,cumsum)
+  
+  if(!fixedflag){
+    findex <- which(!is.na(unlist(fixed)))
+    eindex <- 1:(na+nb)
+    eindex <- eindex[!eindex %in% findex]
+    temp <- rep(0,na+nb); temp2 <- matrix(0,nrow=na+nb,ncol=na+nb)
+
+    temp[eindex] <- coef; temp[findex] <- fixedpars; coef <- temp
+    temp2[eindex,eindex] <- vcov; vcov <- temp2
+  } 
+  
+  if(na==0){
+    A <- 1
+  } else {
+    A <- c(1,coef[1:na])
+  }
+  model <- idpoly(A = A,B = coef[na+1:nb],
+               ioDelay = nk,Ts=deltat(x),noiseVar = sqrt(sigma2),
+               intNoise=intNoise,unit=x$unit)
+  
+  estpoly(sys = model,stats=list(vcov = vcov, sigma = sqrt(sigma2),
+              df = df),fitted.values=fit,residuals=eps[1:N,,drop=F],
+          call=match.call(),input=u)
+}
+
+#' Estimate ARMAX Models
+#' 
+#' Fit an ARMAX model of the specified order given the input-output data 
+#' 
+#' @param x an object of class \code{idframe}
+#' @param order Specification of the orders: the four integer components 
+#' (na,nb,nc,nk) are the order of polynolnomial A, order of polynomial B 
+#' + 1, order of the polynomial C,and the input-output delay respectively
+#' @param init_sys Linear polynomial model that configures the initial parameterization.
+#' Must be an ARMAX model. Overrules the \code{order} argument
+#' @param intNoise Logical variable indicating whether to add integrators in
+#' the noise channel (Default=\code{FALSE})
+#' @param options Estimation Options, setup using \code{\link{optimOptions}}
+#' 
+#' @details
+#' SISO ARMAX models are of the form 
+#' \deqn{
+#'    y[k] + a_1 y[k-1] + \ldots + a_{na} y[k-na] = b_{nk} u[k-nk] + 
+#'    \ldots + b_{nk+nb} u[k-nk-nb] + c_{1} e[k-1] + \ldots c_{nc} e[k-nc]
+#'    + e[k] 
+#' }
+#' The function estimates the coefficients using non-linear least squares 
+#' (Levenberg-Marquardt Algorithm)
+#' \cr
+#' The data is expected to have no offsets or trends. They can be removed 
+#' using the \code{\link{detrend}} function. 
+#' 
+#' @return
+#'  An object of class \code{estpoly} containing the following elements:
+#'  \item{sys}{an \code{idpoly} object containing the 
+#'    fitted ARMAX coefficients}
+#'  \item{fitted.values}{the predicted response}
+#'  \item{residuals}{the residuals}
+#'  \item{input}{the input data used}
+#'  \item{call}{the matched call}
+#'  \item{stats}{A list containing the following fields: \cr
+#'      \code{vcov} - the covariance matrix of the fitted coefficients \cr
+#'      \code{sigma} - the standard deviation of the innovations}
+#'  \item{options}{Option set used for estimation. If no 
+#'    custom options were configured, this is a set of default options}
+#'  \item{termination}{Termination conditions for the iterative
+#'     search used for prediction error minimization:
+#'      \code{WhyStop} - Reason for termination \cr
+#'      \code{iter} - Number of Iterations \cr
+#'      \code{iter} - Number of Function Evaluations }
+#' 
+#' 
+#' @references
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+#' Theory and Practice}, CRC Press, Boca Raton. Sections 14.4.1, 21.6.2
+#' 
+#' @examples
+#' data(armaxsim)
+#' z <- dataSlice(armaxsim,end=1533) # training set
+#' mod_armax <- armax(z,c(1,2,1,2))
+#' mod_armax
+#' 
+#' @export
+armax <- function(x,order=c(0,1,1,0),init_sys=NULL,intNoise=FALSE,
+                  options=optimOptions()){
+  y <- outputData(x); u <- inputData(x)
+  if(intNoise){
+    y <- apply(y,2,diff)
+    u <- apply(u,2,diff)
+  }
+  N <- dim(y)[1]
+  if(!is.null(init_sys)){
+    checkInitSys(init_sys)
+    
+    # Extract orders from initial guess
+    na <- length(init_sys$A) -1;nb <- length(init_sys$B); 
+    nc <- length(init_sys$C) -1;nk <- init_sys$ioDelay
+    order <- c(na,nb,nc,nk)
+    
+    # Initial guess
+    theta0 <- matrix(params(init_sys))
+    ivs <- matrix(predict(init_sys,x))
+    e_init <- y-ivs
+  } else{
+    na <- order[1];nb <- order[2]; nc <- order[3]; nk <- order[4]
+    
+    if(nc<1) 
+      stop("Error: Not an ARMAX model")
+    
+    # Initial Parameter Estimates
+    mod_arx <- iv4(x,c(na,nb,nk)) # fitting ARX model
+    eps_init <- matrix(resid(mod_arx))
+    mod_ma <- arima(eps_init,order=c(0,0,nc),include.mean = F)
+    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))
+  }
+  
+  nb1 <- nb+nk-1 ; n <- max(na,nb1,nc); df <- N - na - nb - nc
+  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))
+  fit <- matrix(y-e)
+  if(intNoise) fit <- apply(fit,2,cumsum)
+  
+  model <- idpoly(A = c(1,theta[1:na]),B = theta[na+1:nb],
+                  C = c(1,theta[na+nb+1:nc]),ioDelay = nk,Ts=deltat(x),
+                  noiseVar = l$sigma,intNoise=intNoise,unit=x$unit)
+  
+  estpoly(sys = model,stats=list(vcov = l$vcov, sigma = l$sigma),
+          fitted.values=fit,residuals=e,call=match.call(),input=u,
+          options = options,termination = l$termination)
+}
+
+#' Estimate Output-Error Models
+#' 
+#' Fit an output-error model of the specified order given the input-output data 
+#' 
+#' @param x an object of class \code{idframe}
+#' @param order Specification of the orders: the four integer components 
+#' (nb,nf,nk) are order of polynomial B + 1, order of the polynomial F,
+#' and the input-output delay respectively
+#' @param init_sys Linear polynomial model that configures the initial parameterization.
+#' Must be an OE model. Overrules the \code{order} argument
+#' @param options Estimation Options, setup using 
+#' \code{\link{optimOptions}}
+#' 
+#' @details
+#' SISO OE models are of the form 
+#' \deqn{
+#'    y[k] + f_1 y[k-1] + \ldots + f_{nf} y[k-nf] = b_{nk} u[k-nk] + 
+#'    \ldots + b_{nk+nb} u[k-nk-nb] + f_{1} e[k-1] + \ldots f_{nf} e[k-nf]
+#'    + e[k] 
+#' }
+#' The function estimates the coefficients using non-linear least squares 
+#' (Levenberg-Marquardt Algorithm)
+#' \cr
+#' The data is expected to have no offsets or trends. They can be removed 
+#' using the \code{\link{detrend}} function. 
+#' 
+#' @return
+#'  An object of class \code{estpoly} containing the following elements:
+#'  \item{sys}{an \code{idpoly} object containing the 
+#'    fitted OE coefficients}
+#'  \item{fitted.values}{the predicted response}
+#'  \item{residuals}{the residuals}
+#'  \item{input}{the input data used}
+#'  \item{call}{the matched call}
+#'  \item{stats}{A list containing the following fields: \cr
+#'      \code{vcov} - the covariance matrix of the fitted coefficients \cr
+#'      \code{sigma} - the standard deviation of the innovations}
+#'  \item{options}{Option set used for estimation. If no 
+#'    custom options were configured, this is a set of default options}
+#'  \item{termination}{Termination conditions for the iterative
+#'     search used for prediction error minimization:
+#'      \code{WhyStop} - Reason for termination \cr
+#'      \code{iter} - Number of Iterations \cr
+#'      \code{iter} - Number of Function Evaluations }
+#' 
+#' @references
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+#' Theory and Practice}, CRC Press, Boca Raton. Sections 14.4.1, 17.5.2, 
+#' 21.6.3
+#' 
+#' @examples
+#' data(oesim)
+#' z <- dataSlice(oesim,end=1533) # training set
+#' mod_oe <- oe(z,c(2,1,2))
+#' mod_oe 
+#' plot(mod_oe) # plot the predicted and actual responses
+#' 
+#' @export
+oe <- function(x,order=c(1,1,0),init_sys=NULL,options=optimOptions()){
+  y <- outputData(x); u <- inputData(x); N <- dim(y)[1]
+  if(!is.null(init_sys)){
+    checkInitSys(init_sys)
+    
+    # Extract orders from initial guess
+    nb <- length(init_sys$B); nf <- length(init_sys$F1) -1
+    nk <- init_sys$ioDelay;order <- c(nb,nf,nk)
+    
+    # Initial guess
+    theta0 <- matrix(params(init_sys))
+    ivs <- matrix(predict(init_sys,x))
+    e_init <- y-ivs
+  } else{
+    nb <- order[1];nf <- order[2]; nk <- order[3];
+    nb1 <- nb+nk-1 ; n <- max(nb1,nf);
+    
+    if(nf<1) 
+      stop("Not an OE model")
+    
+    # Initial Model
+    mod_arx <- iv4(x,c(nf,nb,nk)) # fitting ARX model
+    wk <- resid(mod_arx)
+    e_init <- as.numeric(stats::filter(wk,filter=-mod_arx$sys$A[-1],
+                                       method = "recursive"))
+    ivs <- y-e_init
+    theta0 <- matrix(c(mod_arx$sys$B,mod_arx$sys$A[-1]))
+  }
+  nb1 <- nb+nk-1 ; n <- max(nb1,nf);df <- N - nb - nf
+  
+  l <- levbmqdt(y,u,order,ivs,obj=oeGrad,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],F1 = c(1,theta[nb+1:nf]),
+                  ioDelay = nk,Ts=deltat(x),noiseVar = l$sigma,unit=x$unit)
+  
+  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)
+}
+
+#' Estimate Box-Jenkins Models
+#' 
+#' Fit a box-jenkins model of the specified order from input-output data 
+#' 
+#' @param z an \code{idframe} object containing the data
+#' @param order Specification of the orders: the five integer components 
+#' (nb,nc,nd,nf,nk) are order of polynomial B + 1, order of the polynomial C,
+#' order of the polynomial D, order of the polynomial F, and the 
+#' input-output delay respectively
+#' @param init_sys Linear polynomial model that configures the initial parameterization.
+#' Must be a BJ model. Overrules the \code{order} argument
+#' @param options Estimation Options, setup using 
+#' \code{\link{optimOptions}}
+#' 
+#' @details
+#' SISO BJ models are of the form 
+#' \deqn{
+#'    y[k] = \frac{B(q^{-1})}{F(q^{-1})}u[k-nk] + 
+#'    \frac{C(q^{-1})}{D(q^{-1})} e[k]
+#' }
+#' The orders of Box-Jenkins model are defined as follows:
+#' \deqn{
+#'    B(q^{-1}) = b_1 + b_2q^{-1} + \ldots + b_{nb} q^{-nb+1}
+#' }
+#' 
+#' \deqn{
+#'    C(q^{-1}) = 1 + c_1q^{-1} + \ldots + c_{nc} q^{-nc}
+#' }
+#' 
+#' \deqn{
+#'    D(q^{-1}) = 1 + d_1q^{-1} + \ldots + d_{nd} q^{-nd}
+#' }
+#' \deqn{
+#'    F(q^{-1}) = 1 + f_1q^{-1} + \ldots + f_{nf} q^{-nf}
+#' }
+#' 
+#' The function estimates the coefficients using non-linear least squares 
+#' (Levenberg-Marquardt Algorithm)
+#' \cr
+#' The data is expected to have no offsets or trends. They can be removed 
+#' using the \code{\link{detrend}} function. 
+#' 
+#' @return
+#'  An object of class \code{estpoly} containing the following elements:
+#'  \item{sys}{an \code{idpoly} object containing the 
+#'    fitted BJ coefficients}
+#'  \item{fitted.values}{the predicted response}
+#'  \item{residuals}{the residuals}
+#'  \item{input}{the input data used}
+#'  \item{call}{the matched call}
+#'  \item{stats}{A list containing the following fields: \cr
+#'      \code{vcov} - the covariance matrix of the fitted coefficients \cr
+#'      \code{sigma} - the standard deviation of the innovations}
+#'  \item{options}{Option set used for estimation. If no 
+#'    custom options were configured, this is a set of default options}
+#'  \item{termination}{Termination conditions for the iterative
+#'     search used for prediction error minimization:
+#'      \code{WhyStop} - Reason for termination \cr
+#'      \code{iter} - Number of Iterations \cr
+#'      \code{iter} - Number of Function Evaluations }
+#' 
+#' @references
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+#' Theory and Practice}, CRC Press, Boca Raton. Sections 14.4.1, 17.5.2, 
+#' 21.6.3
+#' 
+#' @examples
+#' data(bjsim)
+#' z <- dataSlice(bjsim,end=1500) # training set
+#' mod_bj <- bj(z,c(2,1,1,1,2))
+#' mod_bj 
+#' residplot(mod_bj) # residual plots
+#' 
+#' @export
+bj <- function(z,order=c(1,1,1,1,0),
+               init_sys=NULL,options=optimOptions()){
+  y <- outputData(z); u <- inputData(z); N <- dim(y)[1]
+  if(!is.null(init_sys)){
+    checkInitSys(init_sys)
+    
+    # Extract orders from initial guess
+    nb <- length(init_sys$B); nf <- length(init_sys$F1) -1
+    nc <- length(init_sys$C) -1;nd <- length(init_sys$d) -1
+    nk <- init_sys$ioDelay;order <- c(nb,nc,nd,nf,nk)
+    
+    # Initial guess
+    theta0 <- matrix(params(init_sys))
+    ivs <- matrix(predict(init_sys,z))
+    e_init <- y-ivs
+  } else{
+    nb <- order[1];nc <- order[2]; nd <- order[3];
+    nf <- order[4]; nk <- order[5];
+    
+    if(nc==0 && nd==0){
+      oe(z,c(nb,nf,nk))
+    } else{
+      
+      # Initial Guess
+      mod_oe <- oe(z,c(nb,nf,nk))
+      v <- resid(mod_oe); zeta <- matrix(predict(mod_oe))
+      mod_arma <- arima(v,order=c(nd,0,nc),include.mean = F)
+      C_params <- if(nc==0) NULL else coef(mod_arma)[nd+1:nc]
+      theta0 <- matrix(c(mod_oe$sys$B,C_params,
+                         -coef(mod_arma)[1:nd],mod_oe$sys$F1[-1]))
+      eps <- matrix(resid(mod_arma))
+    }
+  }
+  l <- levbmqdt(y,u,order,zeta,eps,obj=bjGrad,theta0=theta0,N=N,
+                opt=options)
+  theta <- l$params
+  e <- ts(l$residuals,start = start(y),deltat = deltat(y))
+  
+  C_params <- if(nc==0) NULL else theta[nb+1:nc]
+  model <- idpoly(B = theta[1:nb],C=c(1,C_params),
+                  D=c(1,theta[nb+nc+1:nd]),
+                  F1 = c(1,theta[nb+nc+nd+1:nf]),
+                  ioDelay = nk,Ts=deltat(z),noiseVar = l$sigma,unit=z$unit)
+  
+  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
diff --git a/R/idframe.R b/R/idframe.R
new file mode 100644
index 0000000..e619fca
--- /dev/null
+++ b/R/idframe.R
@@ -0,0 +1,280 @@
+#' S3 class for storing input-output data. 
+#'
+#' \code{idframe} is an S3 class for storing and manipulating input-ouput data. It supports discrete time and frequency domain data.
+#'
+#' @param output dataframe/matrix/vector containing the outputs
+#' @param input dataframe/matrix/vector containing the inputs
+#' @param Ts sampling interval (Default: 1)
+#' @param start Time of the first observation
+#' @param end Time of the last observation Optional Argument
+#' @param unit Time unit (Default: "seconds")
+#'
+#' @return an idframe object
+#' 
+#' @seealso  \code{\link{plot.idframe}}, the plot method for idframe objects
+#' 
+#' @examples
+#' 
+#' dataMatrix <- matrix(rnorm(1000),ncol=5) 
+#' data <- idframe(output=dataMatrix[,3:5],input=dataMatrix[,1:2],Ts=1)
+#' 
+#' @export
+idframe <- function(output,input=NULL,Ts = 1,start=0,end=NULL, 
+                    unit = c("seconds","minutes","hours",
+                             "days")[1]){
+  
+  l <- list(output)
+  if(!is.null(input)) l[[2]] <- input
+  l2 <- lapply(l,data.frame)
+  n <- dim(l2[[1]])
+  dims <- sapply(l2,ncol)
+  colnames(l2[[1]]) <- sapply(1:dims[1],
+                              function(x) paste("y",x,sep = ""))
+  if(!is.null(input)){
+    colnames(l2[[2]]) <- sapply(1:dims[2],
+                              function(x) paste("u",x,sep = "")) 
+  }
+  
+  if(!is.null(end)){
+    start <- end - Ts*(n-1)
+  } 
+  
+  l3 <- lapply(l2,ts,start=start,deltat=Ts)
+  
+  if(!is.null(input)) input <- l3[[2]]
+  # Object Constructor
+  dat <- list(output=l3[[1]],input=input,unit=unit)
+  class(dat) <- "idframe"
+  return(dat)
+}
+
+#' Plotting idframe objects
+#' 
+#' Plotting method for objects inherting from class \code{idframe}
+#' 
+#' @param x an \code{idframe} object
+#' @param col line color, to be passed to plot.(Default=\code{"steelblue"})
+#' @param lwd line width, in millimeters(Default=\code{1})
+#' @param main the plot title. (Default = \code{NULL})
+#' @param size text size (Default = \code{12})
+#' @param \ldots additional arguments
+#' 
+#' @examples
+#' data(cstr)
+#' plot(cstr,col="blue")
+#' 
+#' @import ggplot2 reshape2
+#' 
+#' @export
+plot.idframe <- function(x,col="steelblue",lwd=1,main=NULL,size=12,...){
+  loadNamespace("ggplot2")
+    if(nInputSeries(x)==0){
+      data <- outputData(x)
+    } else{
+      data <- cbind(outputData(x),inputData(x))
+      colnames(data) <- c(outputNames(x),inputNames(x))
+    }
+    df <- melt(data.frame(time=as.numeric(time(data)), data), id.vars="time")
+    with(df,ggplot(df, aes(time, value)) + geom_line(size=lwd,color=col) +
+      facet_grid(variable ~ .,scales="free") + theme_bw(size) + 
+      ylab("Amplitude") + ggtitle(main) + 
+      xlab(paste("Time (",x$unit,")",sep="")))
+      #theme(axis.title.x=element_text(size=11)) + ggtitle(main)
+}
+
+#' @export
+summary.idframe <- function(object,...){
+  out_sum <- summary(outputData(object))
+  in_sum <- summary(inputData(object))
+  
+  out <- list(out_sum=out_sum,in_sum=in_sum,Ts=deltat(object),
+              unit=object$unit,nsample = dim(outputData(object))[1])
+  
+  class(out) <- "summary.idframe"
+  return(out)
+}
+
+#' @export
+print.summary.idframe <- function(x,...){
+  cat("\t\t Number of samples:");cat(x$nsample)
+  cat("\nSampling time: ")
+  cat(x$Ts);cat(" ");cat(x$unit)
+  
+  cat("\n\n")
+  cat("Outputs \n")
+  print(x$out_sum)
+  cat("\n")
+  
+  cat("Inputs \n")
+  print(x$in_sum)
+}
+
+#' Sampling times of IO data
+#'  
+#' \code{time} creates the vector of times at which data was sampled. \code{frequency} returns the number of damples per unit time and \code{deltat} the time-interval 
+#' between observations
+#'  
+#' @param x a idframe object, or a univariate or multivariate time-series, or a vector or matrix
+#'  
+#' @aliases frequency deltat
+#'  
+#' @export
+time <- function(x){
+  if(class(x)[1]!="idframe"){
+    stats::time(x)
+  } else{
+    stats::time(x$output)
+  }
+}
+
+#' @export
+frequency <- function(x){
+  if(class(x)[1]!="idframe"){
+    stats::frequency(x)
+  } else{
+    stats::frequency(x$output)
+  }
+}
+
+#' @export
+deltat <- function(x){
+  if(class(x)[1]!="idframe"){
+    stats::deltat(x)
+  } else{
+    stats::deltat(x$output)
+  }
+}
+
+#' S3 class constructor for storing frequency response data
+#' 
+#' @param respData frequency response data. For SISO systems, supply a 
+#' vector of frequency response values. For MIMO systems with Ny 
+#' outputs and Nu inputs, supply an array of size c(Ny,Nu,Nw).
+#' @param freq frequency points of the response
+#' @param Ts sampling time of data
+#' @param spec power spectra and cross spectra of the system 
+#' output disturbances (noise). Supply an array of size (Ny,Ny,Nw)
+#' @param covData response data covariance matrices. Supply an array
+#' of size (Ny,Nu,Nw,2,2). covData[ky,ku,kw,,] is the covariance matrix
+#' of respData[ky,ku,kw]
+#' @param noiseCov power spectra variance. Supply an array of 
+#' size (Ny,Ny,Nw)
+#' 
+#' @return an idfrd object
+#' 
+#' @seealso
+#' \code{\link{plot.idfrd}} for generating bode plots, 
+#' \code{\link{spa}} and \code{\link{etfe}} for estimating the 
+#' frequency response given input/output data
+#' 
+#' @export
+idfrd <- function(respData,freq,Ts,spec=NULL,covData=NULL,
+                  noiseCov=NULL){
+  # For SISO systems
+  if(is.vector(respData)){
+    dim(respData) <- c(1,1,nrow(freq))
+    if(!is.null(spec)) dim(spec) <- c(1,1,nrow(freq))
+  }
+  
+  out <- list(response=respData,freq=freq,Ts=Ts,spec=spec,covData=
+                covData,noiseCov = noiseCov)
+  class(out) <- "idfrd"
+  return(out)
+}
+
+#' Plotting idfrd objects
+#' 
+#' Generates the bode plot of the given frequency response data. It uses the
+#' ggplot2 plotting engine
+#' 
+#' @param x An object of class \code{idframe}
+#' @param col a specification for the line colour (Default : \code{"
+#' steelblue"})
+#' @param lwd the line width, a positive number, defaulting to 1
+#' @param \ldots additional arguments
+#' 
+#' @seealso \code{\link[ggplot2]{ggplot}}
+#' 
+#' @examples
+#' data(frd)
+#' plot(frd)
+#' 
+#' @import ggplot2 reshape2 signal
+#' 
+#' @export
+plot.idfrd <- function(x,col="steelblue",lwd=1,...){
+  loadNamespace("ggplot2")
+  nfreq <- dim(x$freq)[1]
+  mag <- 20*log10(Mod(x$resp))
+  nout <- dim(mag)[1]; nin <- dim(mag)[2]
+  dim(mag) <- c(nin*nout,nfreq)
+  
+  temp <- aperm(Arg(x$resp),c(3,2,1));dim(temp) <- c(nfreq,nin*nout)
+  l <- t(split(temp, rep(1:ncol(temp), each = nrow(temp))))
+  phase <- 180/pi*t(sapply(l,signal::unwrap))
+
+  g <- vector("list",nin*nout)
+  
+  for(i in 1:length(g)){
+    df <- data.frame(Frequency=x$freq[,],Magnitude=mag[i,],
+                     Phase = phase[i,])
+    melt_df <- reshape2::melt(df,id.var="Frequency")
+    yindex <- (i-1)%/%nin + 1;uindex <- i-nin*(yindex-1)
+    subtitle <- paste("From: u",uindex," to y",yindex,sep="")
+    g[[i]] <- with(melt_df,ggplot(melt_df, aes(Frequency, value)) + 
+      geom_line(size=lwd,color=col) + scale_x_log10(expand = c(0.01,0.01)) + 
+      facet_grid(variable ~ .,scales="free_y") +
+      theme_bw(14) + ylab("") + ggtitle(subtitle) +
+      xlab(ifelse(yindex==nout,"Frequency","")) + 
+      theme(axis.title.x=element_text(color = "black",face = "plain"),
+            title=element_text(size=9,color = "black",face="bold")) + 
+      geom_vline(xintercept=max(x$freq),size=1))
+  }
+  
+  multiplot(plotlist=g,layout=matrix(1:length(g),nrow=nout,byrow=T))
+}
+
+# Multiple plot function
+#
+# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
+# - cols:   Number of columns in layout
+# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
+#
+# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
+# then plot 1 will go in the upper left, 2 will go in the upper right, and
+# 3 will go all the way across the bottom.
+#
+multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
+  
+  # Make a list from the ... arguments and plotlist
+  plots <- c(list(...), plotlist)
+  
+  numPlots = length(plots)
+  
+  # If layout is NULL, then use 'cols' to determine layout
+  if (is.null(layout)) {
+    # Make the panel
+    # ncol: Number of columns of plots
+    # nrow: Number of rows needed, calculated from # of cols
+    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
+                     ncol = cols, nrow = ceiling(numPlots/cols))
+  }
+  
+  if (numPlots==1) {
+    print(plots[[1]])
+    
+  } else {
+    # Set up the page
+    grid::grid.newpage()
+    grid::pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
+    
+    # Make each plot, in the correct location
+    for (i in 1:numPlots) {
+      # Get the i,j matrix positions of the regions that contain this subplot
+      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
+      
+      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
+                                      layout.pos.col = matchidx$col))
+    }
+  }
+}
diff --git a/R/idinput.R b/R/idinput.R
new file mode 100644
index 0000000..46b8cad
--- /dev/null
+++ b/R/idinput.R
@@ -0,0 +1,156 @@
+#' function to generate input singals (rgs/rbs/prbs/sine)
+#' 
+#' \code{idinput} is a function for generating input signals (rgs/rbs/prbs/sine) for identification purposes
+#' 
+#' @param n integer length of the input singal to be generated
+#' @param type the type of input signal to be generated. 
+#' 'rgs' - generates random gaussian signal
+#' 'rbs' - generates random binary signal
+#' 'prbs' - generates pseudorandom binary signal
+#' 'sine' - generates a signal that is a sum of sinusoids
+#' 
+#' Default value is type='rgs'
+#' @param band determines the frequency content of the signal. 
+#' For type='rbs'/'sine'/,  band = [wlow,whigh]
+#' which specifies the lower and the upper bound of the passband frequencies(expressed as fractions of Nyquist frequency). Default is c(0,1)
+#' For type='prbs', band=[0,B]
+#' where B is such that the singal is constant over 1/B (clock period). Default is c(0,1)
+#' @param levels row vector defining the input level. It is of the form 
+#' levels=c(minu, maxu)
+#' For 'rbs','prbs', 'sine', the generated signal always between minu and maxu.
+#' For 'rgs', minu=mean value of signal minus one standard deviation and maxu=mean value of signal plus one standard deviation
+#' 
+#' Default value is levels=c(-1,1)
+#' 
+#' @export
+idinput<-function(n,type='rgs',band=c(0,1),levels=c(-1,1)){
+  if(type=="rbs"){
+    rbs(n,band,levels)
+  } else if(type=="rgs"){
+    rgs(n,band,levels)
+  } else if(type=="sine"){
+    multisine(n,1,band,levels)
+  }else if(type=="prbs"){
+    idin.prbs(n,band,levels)
+  }
+}
+
+rgs <- function(n,band,levels){
+  u <- butter_filt(rnorm(n),band)
+  mu<-(levels[1]+levels[2])/2
+  sigma<-(levels[2]-levels[1])/2
+  u*sigma+mu
+}
+
+rbs <- function(n,band,levels){
+  u <- butter_filt(rnorm(n),band)
+  sapply(u,function(x) if(x>0) levels[2] else levels[1])
+}
+
+butter_filt <- function(x,band){
+  filt <- T; type <- "pass"
+  if(band[1]<=2e-3){
+    if(band[2]==1){
+      filt <- F
+    } else{
+      type <- "low"
+    }
+  } else{
+    if(band[2]==1){
+      type <- "high"  
+    }
+  }
+  if(filt==T){
+    if(type=="low"){
+      bf <- signal::butter(8,band[2],type,"z")
+    } else if(type=="pass"){
+      bf <- signal::butter(8,band,type,"z")
+    }else{
+      bf <- signal::butter(8,band[1],type,"z")
+    }
+    x <- as.numeric(signal::filter(bf,x))
+  }
+  return(matrix(x,ncol=1))
+}
+
+multisine <- function(N,nin=1,band,levels){
+  sinedata <- list(nSin=10,nTrial=10,gridSkip=1)
+  freq <- 2*pi*seq(1,floor(N/2),by=sinedata$gridSkip)/N
+  band <- band*pi
+  freq <- freq[freq>=band[1] & freq<=band[2]]
+  nl <- length(freq)
+  freqInd <- seq(from=1,to=nl,length.out = nin*sinedata$nSin)
+  
+  # Begin Trials
+  trials <- function(x){
+    freqs <- freq[freqInd[seq(from=x,to=nin*sinedata$nSin,by=nin)]]
+    for(k in 1:sinedata$nTrial){
+      utrial <- rowSums(cos(sapply(freqs,
+                                   function(x) (1:N-1)*x+2*pi*rnorm(1))))
+      if(k==1) u <- utrial; temp <- diff(range(utrial))
+      if(diff(range(utrial))< temp) {
+        u <- utrial; temp <- diff(range(utrial))
+      }
+    }
+    clevel <- max(abs(u))
+    diff(levels)/2/clevel*(u+clevel)+levels[1]
+  }
+  u <- sapply(1:nin,trials) 
+  return(u)
+}
+
+#' @import bitops signal
+idin.prbs<-function(n,band=c(0,1),levels=c(0,1)){
+  u <- numeric(0)
+  for(i in 1:18){
+    if(n / (2^i) < 1){
+      u <- idin.prbs12(i,band,levels)
+      break
+    }
+  }
+  return(u[1:n])
+}
+
+idin.prbs12 <- function(N,band=c(0,1),levels=c(0,1)){
+  first=ceiling(abs(rnorm(1)*10))  #some non-zero initial state 
+  x= first
+  v = vector()  
+  n=2^N-1
+  i=1
+  clock=floor(1/band[2])
+  k=1
+  M=rbind(c(0,0,0,0),c(1,2,0,0),c(1,3,0,0),c(1,4,0,0),c(2,5,0,0),c(1,6,0,0),
+          c(1,7,0,0),c(1,2,7,8),c(4,9,0,0),c(3,10,0,0),c(9,11,0,0),
+          c(6,8,11,12),c(9,10,12,13),c(4,8,13,14),c(14,15,0,0),c(4,13,15,16),
+          c(14,17,0,0),c(11,18,0,0))
+  repeat{
+    a=M[N,1]
+    b=M[N,2]
+    c=M[N,3]
+    d=M[N,4]
+    four=c(8,12,13,14,16)
+    if(N %in% four){
+      e=bitwXor(bitwShiftR(x,N-a),bitwShiftR(x,N-b))
+      f=bitwXor(bitwShiftR(x,N-c),bitwShiftR(x,N-d))
+      newbit=bitwAnd(bitwXor(e,f),1)
+    }else{
+      newbit=bitwAnd(bitwXor(bitwShiftR(x,N-a),bitwShiftR(x,N-b)),1)
+    }
+    if(k>=clock || i==1){
+      v[i]=newbit
+      i=i+1
+      k=1
+    }else{
+      v[i]=v[i-1]
+      i=i+1
+      k=k+1
+    }
+    x=bitwOr(bitwShiftR(x,1),bitwShiftL(newbit,N-1))
+    if(x==first){
+      break
+    }
+  }
+  
+  v=sapply(v, function(x){if (x==0) levels[1] else levels[2]})
+  return(v)
+}
\ No newline at end of file
diff --git a/R/imports.R b/R/imports.R
new file mode 100644
index 0000000..e994228
--- /dev/null
+++ b/R/imports.R
@@ -0,0 +1,5 @@
+#' @importFrom stats predict resid arima ts start ar acf ccf coef fitted lm mvfft rnorm window
+#' @importFrom graphics par plot
+#' @importFrom utils read.table tail
+#' @importFrom grid viewport pushViewport grid.layout grid.newpage
+NULL
\ No newline at end of file
diff --git a/R/ioNamesData.R b/R/ioNamesData.R
new file mode 100644
index 0000000..f28f7ae
--- /dev/null
+++ b/R/ioNamesData.R
@@ -0,0 +1,129 @@
+#' Number of series in input or output
+#' 
+#' Number of series in input or output in a idframe object
+#' 
+#' 
+#' @param data \code{idframe} object
+#' 
+#' @aliases nOutputSeries
+#' 
+#' @export
+nInputSeries <- function(data) {
+  ifelse(is.null(data$input),0,ncol(data$input))
+}
+
+#' @export
+nOutputSeries <- function(data) ncol(data$output)
+
+#' Output or Input-data 
+#' 
+#' Extract output-data or input-data in idframe objects
+#' 
+#' @param x \code{idframe} object
+#' @param series the indices to extract
+#' 
+#' @aliases inputData.idframe outputData outputData.idframe
+#' 
+#' @export
+inputData <- function(x,series) UseMethod("inputData")
+
+#' @export
+inputData.default <- function(x,series=NULL){
+  print("Not defined for this datatype")
+}
+
+#' @import tframe
+#' @export
+inputData.idframe <- function(x,series=seq(nInputSeries(x))){
+  if (is.null(x$input))
+    NULL
+  else selectSeries(x$input,series=series)
+}
+
+"inputData<-" <- function(x,value) UseMethod("inputData<-")
+
+"inputData<-.idframe" <- function(x,value){
+  x$input <- value
+  x
+}
+
+
+#' @export
+inputNames <- function(x) UseMethod("inputNames")
+
+#' @export
+inputNames.default <- function(x){
+  print("Not defined for this datatype")
+}
+
+#' @import tframe
+#' @export
+inputNames.idframe <- function(x){
+  seriesNames(inputData(x))
+}
+
+#' Extract or set series' names
+#' 
+#' Extract or set names of series in input or output
+#' 
+#' @param x \code{idframe} object
+#' @param value vector of strings
+#' 
+#' @aliases inputNames inputNames<-.idframe outputNames outputNames outputNames<- outputNames<-.idframe
+#' 
+#' @export
+"inputNames<-" <- function(x,value) UseMethod("inputNames<-")
+
+#' @import tframe
+#' @export
+"inputNames<-.idframe" <- function(x,value){
+  seriesNames(inputData(x)) <- value
+  x
+}
+
+#' @export
+outputData <- function(x,series) UseMethod("outputData")
+
+#' @export
+outputData.default <- function(x,series=NULL){
+  print("Not defined for this datatype")
+}
+
+#' @import tframe
+#' @export
+outputData.idframe <- function(x,series=seq(nOutputSeries(x))){
+  if (is.null(x$output))
+    NULL
+  else selectSeries(x$output,series=series)
+}
+
+"outputData<-" <- function(x,value) UseMethod("outputData<-")
+
+"outputData<-.idframe" <- function(x,value){
+  x$output <- value
+  x
+}
+
+#' @export
+outputNames <- function(x) UseMethod("outputNames")
+
+#' @export
+outputNames.default <- function(x){
+  print("Not defined for this datatype")
+}
+
+#' @import tframe
+#' @export
+outputNames.idframe <- function(x){
+  seriesNames(outputData(x))
+}
+
+#' @export 
+"outputNames<-" <- function(x,value) UseMethod("outputNames<-")
+
+#' @import tframe
+#' @export
+"outputNames<-.idframe" <- function(x,value){
+  seriesNames(outputData(x)) <- value
+  x
+}
diff --git a/R/iv.R b/R/iv.R
new file mode 100644
index 0000000..73b7cc8
--- /dev/null
+++ b/R/iv.R
@@ -0,0 +1,157 @@
+#' ARX model estimation using instrumental variable method
+#' 
+#' Estimates an ARX model of the specified order from input-output data using
+#' the instrument variable method. If arbitrary instruments are not supplied 
+#' by the user, the instruments are generated using the arx routine
+#' 
+#' @param z an idframe object containing the data
+#' @param order Specification of the orders: the three integer components 
+#' (na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) 
+#' and the input-output delay
+#' @param x instrument variable matrix. x must be of the same size as the output 
+#' data. (Default: \code{NULL})
+#' 
+#' @return
+#' An object of class \code{estpoly} containing the following elements:
+#'  \item{sys}{an \code{idpoly} object containing the 
+#'    fitted ARX coefficients}
+#'  \item{fitted.values}{the predicted response}
+#'  \item{residuals}{the residuals}
+#'  \item{input}{the input data used}
+#'  \item{call}{the matched call}
+#'  \item{stats}{A list containing the following fields: \cr
+#'      \code{vcov} - the covariance matrix of the fitted coefficients \cr
+#'      \code{sigma} - the standard deviation of the innovations\cr
+#'      \code{df} - the residual degrees of freedom}
+#' 
+#' @references
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+#' Theory and Practice}, CRC Press, Boca Raton. Sections 21.7.1, 21.7.2
+#' 
+#' Lennart Ljung (1999), \emph{System Identification: Theory for the User}, 
+#' 2nd Edition, Prentice Hall, New York. Section 7.6
+#' 
+#' @examples
+#' data(arxsim)
+#' mod_iv <- iv(arxsim,c(2,1,1))
+#' 
+#' @seealso \code{\link{arx}}, \code{\link{iv4}}
+#' 
+#' @export
+iv <- function(z,order=c(0,1,0),x=NULL){
+  if(is.null(x)){
+    # Initial Guess using ARX
+    mod_arx <- arx(z,order)
+    x <- matrix(sim(mod_arx$sys,inputData(z)))
+  }
+  
+  ivcompute(z,x,order)
+}
+
+#' @import signal
+ivcompute <- function(z,x,order,filt=NULL){
+  y <- outputData(z); u <- inputData(z); N <- dim(y)[1]
+  na <- order[1];nb <- order[2]; nk <- order[3]
+  
+  nb1 <- nb+nk-1 ; n <- max(na,nb1); df <- N-na-nb
+  yout <- apply(y,2,padZeros,n=n);
+  xout <- apply(x,2,padZeros,n=n);
+  uout <- apply(u,2,padZeros,n=n);
+  
+  # Regressors
+  reg <- function(i) {
+    if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1
+    c(-yout[i-1:na,,drop=T],uout[v,,drop=T])
+  }
+  phi <- t(sapply(n+1:(N+n),reg))
+  Y <- yout[n+1:(N+n),,drop=F]
+  
+  # Generating IVs
+  ivx <- function(i) {
+    if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1
+    c(-xout[i-1:na,,drop=T],uout[v,,drop=T])
+  }
+  psi <- t(sapply(n+1:(N+n),ivx))
+  
+  l <- list(phi,psi,Y)
+  if(!is.null(filt)){
+    appfilt <- function(x) apply(x,2,signal::filter,filt=filt)
+    l <- lapply(l, appfilt)
+  }
+  phif <- l[[1]]; psif <- l[[2]]; Yf <- l[[3]]
+  
+  # Estimator
+  lhs <- t(psif)%*%phif; lhsinv <- solve(lhs)
+  theta <- lhsinv%*%t(psif)%*%Yf
+  
+  # Residuals
+  ypred <- (phi%*%theta)[1:N,,drop=F] 
+  e <- y-ypred
+  sigma2 <- norm(e,"2")^2/df
+  vcov <- sigma2*solve(t(phi)%*%phi)
+  
+  model <- idpoly(A = c(1,theta[1:na]),B = theta[na+1:nb],ioDelay = nk,
+                  Ts=deltat(z),noiseVar = sqrt(sigma2),unit=unit)
+  
+  estpoly(sys = model,
+          stats=list(vcov = vcov, sigma = sqrt(sigma2),df = df),
+          fitted.values=ypred,residuals=e,call=match.call(),input=u)
+}
+
+#' ARX model estimation using four-stage instrumental variable method
+#' 
+#' Estimates an ARX model of the specified order from input-output data using
+#' the instrument variable method. The estimation algorithm is insensitive to 
+#' the color of the noise term.
+#' 
+#' @param z an idframe object containing the data
+#' @param order Specification of the orders: the three integer components 
+#' (na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) 
+#' and the input-output delay
+#' 
+#' @details 
+#' Estimation is performed in 4 stages. The first stage uses the arx function. The resulting model generates the 
+#' instruments for a second-stage IV estimate. The residuals obtained from this model are modeled using a sufficently
+#' high-order AR model. At the fourth stage, the input-output data is filtered through this AR model and then subjected 
+#' to the IV function with the same instrument filters as in the second stage.
+#' 
+#' @return
+#' An object of class \code{estpoly} containing the following elements:
+#'  \item{sys}{an \code{idpoly} object containing the 
+#'    fitted ARX coefficients}
+#'  \item{fitted.values}{the predicted response}
+#'  \item{residuals}{the residuals}
+#'  \item{input}{the input data used}
+#'  \item{call}{the matched call}
+#'  \item{stats}{A list containing the following fields: \cr
+#'      \code{vcov} - the covariance matrix of the fitted coefficients \cr
+#'      \code{sigma} - the standard deviation of the innovations\cr
+#'      \code{df} - the residual degrees of freedom}
+#' 
+#' @references
+#' Lennart Ljung (1999), \emph{System Identification: Theory for the User}, 
+#' 2nd Edition, Prentice Hall, New York. Section 15.3
+#' 
+#' @examples
+#' mod_dgp <- idpoly(A=c(1,-0.5),B=c(0.6,-.2),C=c(1,0.6),ioDelay = 2,noiseVar = 0.1)
+#' u <- idinput(400,"prbs")
+#' y <- sim(mod_dgp,u,addNoise=TRUE)
+#' z <- idframe(y,u)
+#' mod_iv4 <- iv4(z,c(1,2,2))
+#' 
+#' @seealso \code{\link{arx}}, \code{\link{iv4}}
+#' @export
+iv4 <- function(z,order=c(0,1,0)){
+  na <- order[1]; nb <- order[2]
+  # Steps 1-2
+  mod_iv <- iv(z,order)
+  
+  # Step 3
+  w <- resid(mod_iv)
+  mod_ar <- ar(w,aic = F,order.max =na+nb)
+  Lhat <- signal::Arma(b=c(1,-mod_ar$ar),a=1)
+  
+  # Step 4
+  x2 <- matrix(sim(mod_iv$sys,inputData(z)))
+  ivcompute(z,x2,order,Lhat)
+}
\ No newline at end of file
diff --git a/R/nonparam.R b/R/nonparam.R
new file mode 100644
index 0000000..8e8e412
--- /dev/null
+++ b/R/nonparam.R
@@ -0,0 +1,271 @@
+#' Estimate Impulse Response Coefficients
+#' 
+#' \code{impulseest} is used to estimate impulse response coefficients from 
+#' the data
+#' 
+#' @param x an object of class \code{idframe}
+#' @param M Order of the FIR Model (Default:\code{30})
+#' @param K Transport delay in the estimated impulse response 
+#' (Default:NULL)
+#' @param regul Parameter indicating whether regularization should be 
+#' used. (Default:\code{FALSE})
+#' @param lambda The value of the regularization parameter. Valid only if
+#' \code{regul=TRUE}. (Default:\code{1})
+#' 
+#' @details
+#' The IR Coefficients are estimated using linear least squares. Future 
+#' Versions will provide support for multivariate data.
+#' 
+#' @references
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+#' Theory and Practice}, CRC Press, Boca Raton. Sections 17.4.11 and 20.2
+#' 
+#' @seealso \code{\link{step}}
+#' 
+#' @examples
+#' uk <- rnorm(1000,1)
+#' yk <- filter (uk,c(0.9,-0.4),method="recursive") + rnorm(1000,1)
+#' data <- idframe(output=data.frame(yk),input=data.frame(uk))
+#' fit <- impulseest(data)
+#' impulseplot(fit)
+#' 
+#' @export
+impulseest <- function(x,M=30,K=NULL,regul=F,lambda=1){
+  
+  N <- dim(x$output)[1]
+  if(is.null(K))
+    K <- rep(0,nInputSeries(x)*nOutputSeries(x))
+  
+  out <- rep(list(0),length(K))
+  
+  for(i in seq(nOutputSeries(x))){
+    for(j in seq(nInputSeries(x))){
+      index <- (i-1)*nInputSeries(x)+j 
+      out[[index]] <- impulsechannel(outputData(x)[,i,drop=F],
+                                     inputData(x)[,j,drop=F],N,M,
+                                     K[index],regul,lambda)
+    }
+  }
+  out$ninputs <- nInputSeries(x)
+  out$noutputs <- nOutputSeries(x)
+  class(out) <- "impulseest"
+  return(out)
+}
+
+impulsechannel <- function(y,u,N,M,K=0,regul=F,lambda=1){
+  ind <- (M+K+1):N 
+  z_reg <- function(i) u[(i-K):(i-M-K),]
+  Z <- t(sapply(ind,z_reg))
+  Y <- y[ind,]
+  
+  # Dealing with Regularization
+  if(regul==F){
+    # Fit Linear Model and find standard errors
+    fit <- lm(Y~Z-1)
+    coefficients <- coef(fit); residuals <- resid(fit)
+  } else{
+    inner <- t(Z)%*%Z + lambda*diag(dim(Z)[2])
+    pinv <- solve(inner)%*% t(Z)
+    coefficients <- pinv%*%Y
+    residuals <- Y - Z%*%coefficients
+  }
+  df <- nrow(Z)-ncol(Z);sigma2 <- sum(residuals^2)/df
+  vcov <- sigma2 * solve(t(Z)%*%Z)
+  se <- sqrt(diag(vcov))
+  
+  out <- list(coefficients=coefficients,residuals=residuals,lags=K:(M+K),
+              x=colnames(u),y=colnames(y),se = se)
+  out
+}
+
+#' Impulse Response Plots
+#' 
+#' Plots the estimated IR coefficients along with the significance limits
+#' at each lag. 
+#' 
+#' @param model an object of class \code{impulseest}
+#' @param sd standard deviation of the confidence region (Default: \code{2})
+#' 
+#' @seealso \code{\link{impulseest}},\code{\link{step}}
+#' @import ggplot2
+#' 
+#' @export
+impulseplot <- function(model,sd=2){
+  plotseq <- seq(model$noutputs*model$ninputs)
+  g <- vector("list",model$nin*model$nout)
+  
+  for(i in plotseq){
+    z <- model[[i]]
+    lim <- z$se*sd
+    yindex <- (i-1)%/%model$nin + 1;uindex <- i-model$nin*(yindex-1)
+    df <- data.frame(x=z$lags,y=coef(z),lim=lim)
+    g[[i]] <- with(df,ggplot(df,aes(x,y))+ 
+      geom_segment(aes(xend=x,yend=0))+geom_hline(yintercept = 0) + 
+      geom_point(size=2) + ggtitle(paste("From",z$x,"to",z$y))+
+      geom_line(aes(y=lim),linetype="dashed",colour="steelblue") +
+      geom_line(aes(y=-lim),linetype="dashed",colour="steelblue") +
+      ggplot2::theme_bw(14) + ylab(ifelse(uindex==1,"IR Coefficients","")) +
+      xlab(ifelse(yindex==model$nout,"Lags","")) + 
+      theme(axis.title=element_text(size=12,color = "black",face = "plain"),
+            title=element_text(size=9,color = "black",face="bold")) +
+      scale_x_continuous(expand = c(0.01,0.01)))
+  }
+  multiplot(plotlist=g,layout=plotseq)
+}
+
+
+#' Step Response Plots
+#' 
+#' Plots the step response of a system, given the IR model
+#' 
+#' @param model an object of class \code{impulseest}
+#' 
+#' @seealso \code{\link{impulseest}}
+#' 
+#' @examples
+#' uk <- rnorm(1000,1)
+#' yk <- filter (uk,c(0.9,-0.4),method="recursive") + rnorm(1000,1)
+#' data <- idframe(output=data.frame(yk),input=data.frame(uk))
+#' fit <- impulseest(data)
+#' step(fit) 
+#' 
+#' @import ggplot2
+#' @export 
+step <- function(model){
+  plotseq <- seq(model$noutputs*model$ninputs)
+  g <- vector("list",model$nin*model$nout)
+  
+  for(i in plotseq){
+    z <- model[[i]]
+    stepResp <- cumsum(coef(z))
+    yindex <- (i-1)%/%model$nin + 1;uindex <- i-model$nin*(yindex-1)
+    df <- data.frame(x=z$lags,y=stepResp)
+    g[[i]] <- with(df,ggplot(df,aes(x,y))+ 
+      geom_step() + ggtitle(paste("From",z$x,"to",z$y)) +
+      theme_bw(14) + ylab(ifelse(uindex==1,"Step Response","")) +
+      xlab(ifelse(yindex==model$nout,"Lags","")) + 
+      theme(axis.title=element_text(size=12,color = "black",face = "plain"),
+            title=element_text(size=9,,color = "black",face="bold")))
+  }
+  multiplot(plotlist=g,layout=plotseq)
+}
+
+#' Estimate frequency response 
+#' 
+#' Estimates frequency response and noise spectrum from data with 
+#' fixed resolution using spectral analysis
+#' 
+#' @param x an \code{idframe} object
+#' @param winsize lag size of the Hanning window (Default: \code{min
+#' (length(x)/10,30)})
+#' @param freq frequency points at which the response is evaluated
+#' (Default: \code{seq(1,128)/128*pi/Ts})
+#' 
+#' @return
+#' an \code{idfrd} object containing the estimated frequency response 
+#' and the noise spectrum
+#' 
+#' @references
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+#' Theory and Practice}, CRC Press, Boca Raton. Sections 16.5 and 20.4
+#' 
+#' @examples
+#' data(arxsim)
+#' frf <- spa(arxsim)
+#' 
+#' @import signal
+#' @export
+spa <- function(x,winsize=NULL,freq=NULL){
+  N <- dim(x$out)[1]
+  nout <- nOutputSeries(x); nin <- nInputSeries(x)
+  
+  if(is.null(winsize)) winsize <- min(N/10,30)
+  if(is.null(freq)) freq <- (1:128)/128*pi/deltat(x)
+  M <- winsize
+  
+  Ryu <- mult_ccf(x$out,x$input,lag.max = M)
+  Ruu <- mult_ccf(x$input,x$input,lag.max=M)
+  Ryy <- mult_ccf(x$out,x$out,lag.max = M)
+  
+  cov2spec <- function(omega,R,M){
+    seq1 <- exp(-1i*(-M:M)*omega)
+    sum(R*signal::hanning(2*M+1)*seq1)
+  }
+  
+  G <- array(0,c(nout,nin,length(freq)))
+  spec <- array(0,c(nout,1,length(freq)))
+  for(i in 1:nout){
+    phi_y <- sapply(freq,cov2spec,Ryy[i,i,],M)
+    temp <- phi_y
+    for(j in 1:nin){
+      phi_yu <- sapply(freq,cov2spec,Ryu[i,j,],M)
+      phi_u <- sapply(freq,cov2spec,Ruu[j,j,],M)
+      G[i,j,] <- phi_yu/phi_u
+      temp <- temp - phi_yu*Conj(phi_yu)/phi_u
+    }
+    spec[i,1,] <- temp
+  }
+  out <- idfrd(G,matrix(freq),deltat(x),spec)
+  return(out)
+}
+
+mult_ccf <- function(X,Y=NULL,lag.max=30){
+  N <- dim(X)[1]; nx <- dim(X)[2]
+  ny <- ifelse(is.null(Y),nx,dim(Y)[2])
+  
+  ccvfij <- function(i,j,lag=30) ccf(X[,i],Y[,j],plot=F,lag.max =lag,
+                                     type="covariance")
+  Xindex <- matrix(sapply(1:nx,rep,nx),ncol=1)[,1]
+  temp <- mapply(ccvfij,i=Xindex,j=rep(1:ny,ny),
+                 MoreArgs = list(lag=lag.max))
+  
+  ccfextract <- function(i,l) l[,i]$acf
+  temp2 <- t(sapply(1:(nx*ny),ccfextract,l=temp))
+  dim(temp2) <- c(nx,ny,2*lag.max+1)
+  return(temp2)
+}
+
+#' Estimate empirical transfer function
+#' 
+#' Estimates the emperical transfer function from the data by taking the 
+#' ratio of the fourier transforms of the output and the input variables
+#' 
+#' @param data an object of class \code{idframe}
+#' @param n frequency spacing (Default: \code{128})
+#' 
+#' @return
+#' an \code{idfrd} object containing the estimated frequency response
+#' 
+#' @references
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+#' Theory and Practice}, CRC Press, Boca Raton. Sections 5.3 and 20.4.2
+#' 
+#' @seealso \code{\link[stats]{fft}}
+#' 
+#' @examples
+#' data(arxsim)
+#' frf <- etfe(arxsim)
+#' 
+#' @export
+etfe <- function(data,n=128){
+  y <- data$output
+  u <- data$input
+  N <- dim(data$output)[1]
+  if(N < n){
+    n=N
+  }
+  v=seq(1,N,length.out = n)
+  temp <- cbind(data$output[v,],data$input[v,])
+  tempfft <- mvfft(temp)/dim(temp)[1]
+  G <- comdiv(tempfft[,1],tempfft[,2])
+  resp = G[1:ceiling(length(G)/2)]
+  frequency <- matrix(seq( 1 , ceiling(n/2) ) * pi / floor(n/2) / deltat(data))
+  out <- idfrd(respData = resp,freq=frequency,Ts=data$Ts)
+  return(out)
+}
+comdiv <- function(z1,z2){
+  mag1 <- Mod(z1);mag2 <- Mod(z2)
+  phi1 <- Arg(z1); phi2 <- Arg(z2)
+  
+  complex(modulus=mag1/mag2,argument=signal::unwrap(phi1-phi2))
+}
\ No newline at end of file
diff --git a/R/poly.R b/R/poly.R
new file mode 100644
index 0000000..528a599
--- /dev/null
+++ b/R/poly.R
@@ -0,0 +1,196 @@
+#' Polynomial model with identifiable parameters
+#' 
+#' Creates a polynomial model with identifiable coefficients
+#' 
+#' @param A autoregressive coefficients
+#' @param B,F1 coefficients of the numerator and denominator respectively
+#' of the deterministic model between the input and output 
+#' @param C,D coefficients of the numerator and denominator respectively
+#' of the stochastic model 
+#' @param ioDelay the delay in the input-output channel
+#' @param Ts sampling interval
+#' @param noiseVar variance of the white noise source (Default=\code{1})
+#' @param intNoise Logical variable indicating presence or absence of integrator
+#' in the noise channel (Default=\code{FALSE})
+#' @param unit time unit (Default=\code{"seconds"})
+#' 
+#' @details
+#' Discrete-time polynomials are of the form
+#' \deqn{
+#'  A(q^{-1}) y[k] = \frac{B(q^{-1})}{F1(q^{-1})} u[k] + 
+#'  \frac{C(q^{-1})}{D(q^{-1})} e[k] 
+#' }
+#' 
+#' @examples
+#' # define output-error model
+#' mod_oe <- idpoly(B=c(0.6,-0.2),F1=c(1,-0.5),ioDelay = 2,Ts=0.1,
+#' noiseVar = 0.1)
+#' 
+#' # define box-jenkins model with unit variance
+#' B <- c(0.6,-0.2)
+#' C <- c(1,-0.3)
+#' D <- c(1,1.5,0.7)
+#' F1 <- c(1,-0.5)
+#' mod_bj <- idpoly(1,B,C,D,F1,ioDelay=1)
+#' 
+#' @export
+idpoly <- function(A=1,B=1,C=1,D=1,F1=1,ioDelay=0,Ts=1,
+                   noiseVar=1,intNoise = F,unit = c("seconds","minutes",
+                                       "hours","days")[1]){
+  Bindex <- which.max(B!=0)
+  ioDelay <- ifelse(Bindex-1==0,ioDelay,Bindex-1)
+  B <- B[Bindex:length(B)]
+  out <- list(A= A,B=B,C=C,D=D,F1=F1,ioDelay = ioDelay,Ts=Ts,
+              noiseVar=noiseVar,unit=unit,intNoise = intNoise)
+  out$type <- typecheck(out)
+  class(out) <- "idpoly"
+  return(out)
+}
+
+typecheck <- function(x){
+  y <- lapply(x[1:5],checkUnity)
+  if(y$A){
+    out <- if(y$C && y$D) if(y$F1) "fir" else "oe" else "bj" 
+  } else{
+    if(y$D && y$F1){
+      if(x$intNoise){
+        out <- if(y$C) "ari" else "arima"
+      } else{
+        out <- if(y$C) "ar" else "arma"
+      }
+      if(!y$B) out <- paste(out,"x",sep="")
+    } else{
+      out <- "polynomial"
+    }
+  }
+  out
+}
+
+checkUnity <- function(x){
+  out <- if(length(x)==1 && x==1) TRUE else FALSE
+}
+
+#' @export
+print.idpoly <- function(x,se=NULL,dig=3,...){
+  main <- paste("Discrete-time",toupper(x$type),"model:")
+  if(x$type=="oe" || x$type=="bj"){
+    main <- paste(main,"y[k] = B(z)/F(z) u[k] +")
+    if(x$type=="bj"){
+      polyExp <- if(!checkUnity(x$C)) "C(z)/D(z)" else "1/D(z)"
+      if(x$intNoise==T) polyExp <- paste(polyExp,"1/(1-z^{-1})")
+      main <- paste(main,polyExp,"e[k] \n\n")
+    }
+  } else{
+    main <- paste(main,"A(z)y[k] =")
+    if(!checkUnity(x$B)) main <- paste(main,"B(z)u[k] +")
+    if(checkUnity(x$C)){
+      Cexp <- if(!x$intNoise) "" else "1/(1-z^{-1})"
+    }else{
+      Cexp <- "C(z)"
+      if(x$intNoise) Cexp <- paste(Cexp,"/(1-z^{-1})",sep="")
+    }
+    main <- paste(main,Cexp,"e[k] \n\n")
+  }
+  cat(main)
+  
+  # Printing Standard error sequence
+  j=1
+  print_se <- function(se){
+    if(!is.null(se)){
+      cat(" (+/- ",round(se[j],dig),") ",sep = "")
+      j <<- j+1
+    }
+  }
+  
+  if(length(x$A)>1){
+    cat("A(q^{-1}) = ")
+    for(i in seq_along(x$A)){
+      if(i-1==0){
+        cat(round(x$A[i],dig))
+      } else{
+        if(x$A[i]>0) cat(" + ") else cat("- ")
+        
+        if(!(abs(x$A[i])==1)) cat(abs(round(x$A[i],dig)))
+        print_se(se)
+        cat("q^{-",i-1,"}",sep="")
+      }
+      cat("\t")
+    }
+    cat("\n")
+  }
+  
+  cat("B(q^{-1}) = ")
+  for(i in seq_along(x$B)){
+    if(i+x$ioDelay-1==0){
+      cat(round(x$B[i],dig))
+    } else{
+      
+      if(!((x$ioDelay!=0) && (i==1))){
+        if(x$B[i]>0) cat(" + ") else cat("- ")
+      } else{
+        if(x$B[i]<0) cat("-")
+      }
+      
+      if(!(abs(x$B[i])==1)) cat(abs(round(x$B[i],dig)))
+      print_se(se)
+      cat("q^{-",i+x$ioDelay-1,"}",sep="")
+    }
+    cat("\t")
+  }
+  cat("\n")
+  
+  if(length(x$C)>1){
+    cat("C(q^{-1}) = ")
+    for(i in seq_along(x$C)){
+      if(i-1==0){
+        cat(round(x$C[i],dig))
+      } else{
+        if(x$C[i]>0) cat(" + ") else cat("- ")
+        
+        if(!(abs(x$C[i])==1)) cat(abs(round(x$C[i],dig)))
+        print_se(se)
+        cat("q^{-",i-1,"}",sep="")
+      }
+      cat("\t")
+    }
+    cat("\n")
+  }
+  
+  if(length(x$D)>1){
+    cat("D(q^{-1}) = ")
+    for(i in seq_along(x$D)){
+      if(i-1==0){
+        cat(round(x$D[i],dig))
+      } else{
+        if(x$D[i]>0) cat(" + ") else cat("- ")
+        
+        if(!(abs(x$D[i])==1)) cat(abs(round(x$D[i],dig)))
+        print_se(se)
+        cat("q^{-",i-1,"}",sep="")
+      }
+      cat("\t")
+    }
+    cat("\n")
+  }
+  
+  if(length(x$F1)>1){
+    cat("F(q^{-1}) = ")
+    for(i in seq_along(x$F1)){
+      if(i-1==0){
+        cat(round(x$F1[i],dig))
+      } else{
+        if(x$F1[i]>0) cat(" + ") else cat("- ")
+        
+        if(!(abs(x$F1[i])==1)) cat(abs(round(x$F1[i],dig)))
+        print_se(se)
+        cat("q^{-",i-1,"}",sep="")
+      }
+      cat("\t")
+    }
+  }
+  cat("\n")
+}
+
+params <- function(x){
+  c(x$A[-1],x$B,x$C[-1],x$D[-1],x$F1[-1])
+} 
\ No newline at end of file
diff --git a/R/predict.R b/R/predict.R
new file mode 100644
index 0000000..b56fa8f
--- /dev/null
+++ b/R/predict.R
@@ -0,0 +1,120 @@
+predict.idpoly <- function(object,data,nahead=1){
+  y <- outputData(data); u<- inputData(data)
+  G <- signal::Arma(b=c(rep(0,object$ioDelay),object$B),
+                    a= as.numeric(polynom::polynomial(object$A)*
+                                    polynom::polynomial(object$F1)))
+  det_sys <- as.numeric(signal::filter(G,u))
+  if(object$type=="oe" || nahead==Inf){
+    ypred <- det_sys
+  } else{
+    Hden <- as.numeric(polynom::polynomial(object$A)*polynom::polynomial(object$D))
+    Hinv <- signal::Arma(b=Hden,a=object$C)
+    filtered <- as.numeric(signal::filter(Hinv,as.numeric(y)-det_sys))
+    if(nahead!=1){
+      H <- as.numeric(polynom::polynomial(object$C)*polyinv(Hden,nahead))
+      Hl <- signal::Ma(H[1:nahead])
+      filtered <- as.numeric(signal::filter(Hl,filtered))
+    }
+    ypred <- as.numeric(y) - filtered
+  }
+  
+  matrix(ts(ypred,start=start(data),deltat=deltat(data)))
+}
+
+polyinv <- function(x,k){
+  gamma <- 1/Re(polyroot(x))
+  
+  inverse <- function(y,k){
+    sapply(1:k-1,function(i) y^i)
+  }
+  z <- lapply(lapply(gamma,inverse,k=2),polynom::polynomial)
+  temp = z[[1]]
+  if(length(z)>1){
+    for(i in 2:length(z)){
+      temp = temp*z[[i]]
+    }
+  }
+  temp
+}
+
+#' Predictions of identified model
+#' 
+#' Predicts the output of an identified model (\code{estpoly}) object K steps ahead. 
+#' 
+#' @param object \code{estpoly} object containing the identified model
+#' @param newdata optional dataset to be used for predictions. If not supplied, 
+#' predictions are made on the training set.
+#' @param nahead number of steps ahead at which to predict (Default:1). For infinite-
+#' step ahead predictions or pure simulation, supply \code{Inf}.
+#' @param \ldots other arguments
+#' 
+#' @return 
+#' Time-series containing the predictions
+#' 
+#' @examples 
+#' data(arxsim)
+#' mod1 <- oe(arxsim,c(2,1,1))
+#' Yhat <- predict(mod1,arxsim) #  1-step ahead predictions 
+#' Yhat_2 <- predict(mod1,arxsim,nahead=2) # 2-step ahead predictions
+#' Yhat_inf <- predict(mod1,arxsim,nahead=Inf) # Infinite-step ahead predictions
+#' 
+#' @references 
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: Theory 
+#' and Practice}, CRC Press, Boca Raton. Chapter 18
+#' 
+#' @export
+predict.estpoly <- function(object,newdata=NULL,nahead=1,...){
+  if(is.null(newdata)&& nahead==1){
+    return(matrix(fitted(object)))
+  } else{
+    model <- object$sys
+    if(is.null(newdata)){
+      y <- fitted(object)+resid(object)
+      u <- object$input
+      z <- idframe(y,u,Ts = deltat(y),start=start(y))
+    } else{
+      z <- newdata
+    }
+    predict(model,z,nahead)
+  } 
+}
+
+#' Compare the measured output and the predicted output(s)
+#' 
+#' Plots the output predictions of model(s) superimposed over validation data, 
+#' data, for comparison. 
+#' 
+#' @param data validation data in the form of an \code{idframe} object
+#' @param nahead number of steps ahead at which to predict (Default:1). For infinite-
+#' step ahead predictions, supply \code{Inf}.
+#' @param \ldots models whose predictions are to be compared
+#' 
+#' @examples 
+#' data(arxsim)
+#' mod1 <- arx(arxsim,c(1,2,2))
+#' mod2 <- oe(arxsim,c(2,1,1))
+#' compare(arxsim,nahead=Inf,mod1,mod2)
+#' 
+#' @seealso \code{\link{predict.estpoly}} for obtaining model predictions
+#' @import ggplot2 reshape2
+#' @export
+compare <- function(data,nahead=1,...){
+  loadNamespace("ggplot2")
+  # Input Validation
+  input_list <- as.list(substitute(list(...)))[-1]
+  dots <- list(...)
+  if(is.null(dots)) stop("No model supplied")
+  
+  Y <- sapply(dots,predict,newdata=data,nahead=nahead)
+  nrmse <- sapply(dots,FUN = function(x) fitch(x)$FitPer)
+  temp <- paste(as.character(input_list),paste(round(nrmse,2),"%",sep=""),
+                       sep=": ")
+  df <- data.frame(Time = as.numeric(time(data)),
+                   Actual=as.numeric(outputData(data)[,1]),Y)
+  meltdf <- reshape2::melt(df,id="Time")
+  
+  with(meltdf,ggplot(meltdf,aes(x=Time,y=value,color=variable,group=variable))+geom_line(size=1)+
+    ggtitle(paste(nahead,"step(s) ahead model predictions"))+
+    theme_bw()+ylab(outputNames(data)) + labs(colour="") + 
+    scale_colour_hue(labels=c("Measured",temp),l=50))
+}
\ No newline at end of file
diff --git a/R/preprocess.R b/R/preprocess.R
new file mode 100644
index 0000000..b5ba991
--- /dev/null
+++ b/R/preprocess.R
@@ -0,0 +1,192 @@
+#' Remove offsets and linear trends
+#' 
+#' Removes offsets or trends from data
+#' @param x an object of class \code{idframe}
+#' @param type argument indicating the type of trend to be removed (Default=\code{0})
+#'  \itemize{
+#'    \item type=\code{0}: Subtracts mean value from each signal
+#'    \item type=\code{1}: Subtracts a linear trend (least-squres fit)
+#'    \item type=\code{trInfo} object: Subtracts a trend specified by the object
+#'  }
+#' 
+#' @return 
+#' A list containing two objects: the detrended data and the trend information
+#' 
+#' @details 
+#' 
+#' \code{R} by default doesn't allow return of multiple objects. The \code{\%=\%}
+#' operator and \code{g} function in this package facillitate this behaviour. See 
+#' the examples section for more information.
+#' 
+#' @aliases trInfo
+#' 
+#' @examples
+#' data(cstr)
+#' datatrain <- dataSlice(cstr,end=4500)
+#' datatest <- dataSlice(cstr,4501)
+#' g(Ztrain,tr) %=% detrend(datatrain) # Remove means
+#' g(Ztest) %=% detrend(datatest,tr)
+#' 
+#' @seealso \code{\link[stats]{lm}}
+#' @export
+detrend <- function(x,type=0){ 
+  z <- x 
+  reg <- time(x)
+  if(class(type)=="trInfo"){ # remove custom trend
+    if(nOutputSeries(x)!=0){
+      fit <- sweep(sweep(matrix(rep(reg,nOutputSeries(x)),ncol=nOutputSeries(x)),
+                         2,type$OutputSlope,"*"),2,type$OutputOffset,"+")
+      z$output <- x$output-fit
+    }
+    if(nInputSeries(x)!=0){
+      fit <- sweep(sweep(matrix(rep(reg,nInputSeries(x)),ncol=nInputSeries(x)),
+                         2,type$InputSlope,"*"),2,type$InputOffset,"+")
+      z$input <- x$input-fit
+    }
+    tinfo <- type
+  } else if(type == 0){ # remove means
+    tinfo <- trInfo()
+    if(nOutputSeries(x)!=0){
+      outputData(z) <- apply(outputData(x),2,scale,T,F)
+      tinfo$OutputOffset <- colMeans(x$output)
+      tinfo$OutputSlope <- rep(0,nOutputSeries(x))
+    }
+    
+    if(nInputSeries(x)!=0){
+      inputData(z) <- apply(inputData(x),2,scale,T,F)
+      tinfo$InputOffset <- colMeans(x$input)
+      tinfo$InputSlope <- rep(0,nInputSeries(x))
+    }
+  } else if(type==1){
+    formula <- X ~ reg
+    
+    # Function which performs linear regression across every column
+    multilm <- function(x,formula,time){
+      l <- lapply(as.list(x),function(x) data.frame(X=x,reg=time))
+      trend <- lapply(l,function(x) lm(formula,data=x))
+      trend
+    }
+    
+    tinfo <- trInfo()
+    if(nOutputSeries(x)!=0){
+      output_trend <- multilm(outputData(x),formula,reg)
+      outputData(z) <- ts(sapply(output_trend,resid),start=reg[1],
+                          end=tail(reg,n=1),deltat=deltat(x))
+      out_coefs <- sapply(output_trend,coef)
+      tinfo$OutputOffset <- out_coefs[1,]
+      tinfo$OutputSlope <- out_coefs[2,]
+    }
+    
+    if(nInputSeries(x)!=0){
+      input_trend <- multilm(inputData(x),formula,reg)
+      inputData(z) <- ts(sapply(input_trend,resid),start=reg[1],
+                         end=tail(reg,n=1),deltat=deltat(x))
+      in_coefs <- sapply(input_trend,coef)
+      tinfo$InputOffset <- in_coefs[1,]
+      tinfo$InputSlope <- in_coefs[2,]
+    }
+  } else{
+    stop("Error: Invalid trend type")
+  }
+  list(z,tinfo)
+}
+
+#' @export
+trInfo <- function(InputOffset=numeric(0),OutputOffset=numeric(0),
+                      InputSlope=numeric(0),OutputSlope=numeric(0)){
+  l <- list(InputOffset=InputOffset,OutputOffset=OutputOffset,
+            InputSlope=InputSlope,OutputSlope=OutputSlope)
+  class(l) <- "trInfo"
+  l
+}
+
+#' Replace Missing Data by Interpolation
+#' 
+#' Function for replacing missing values with interpolated ones. This is an
+#' extension of the \code{na.approx} function from the \code{zoo} package.
+#' The missing data is indicated using the value \emph{NA}.
+#' 
+#' @param data an object of class \code{idframe}
+#' @return 
+#' data (an idframe object) with missing data replaced.
+#' 
+#' @seealso \code{\link[zoo]{na.approx}}
+#' 
+#' @examples
+#' data(cstr_mis)
+#' summary(cstr_mis) # finding out the number of NAs
+#' cstr <- misdata(cstr_mis)
+#' 
+#' @importFrom zoo na.approx
+#' @export
+misdata <- function(data){
+  if (!requireNamespace("zoo", quietly = TRUE)) {
+    stop("Package zoo needed for this function to work. Please install it.",
+         call. = FALSE)
+  }
+  
+  f <- function(var,start,end,Ts){
+    time_range <- range(time(var))
+    start <- time_range[1];end <- time_range[2]
+    Ts <- stats::deltat(var)
+    var <- ts(data=var,start=start,end=end,deltat=Ts)
+    out <- zoo::na.approx(var,na.rm=F)
+    return(as.numeric(out))
+  }
+  
+  Z <- data
+  outputData(Z) <- apply(outputData(data),2,f)
+  inputData(Z) <- apply(inputData(data),2,f)
+  Z
+}
+
+
+#' Subset or Resample idframe data
+#' 
+#' \code{dataSlice} is a subsetting method for objects of class \code{idframe}. It 
+#' extracts the subset of the object \code{data} observed between indices \code{start}
+#' and \code{end}.  If a frequency is specified, the series is then re-sampled at the 
+#' new frequency.
+#' 
+#' @param data an object of class \code{idframe}
+#' @param start the start index
+#' @param end the end index
+#' @param freq fraction of the original frequency at which the series
+#' to be sampled.
+#' 
+#' @details
+#' The dataSlice function extends the \code{\link[stats]{window}} 
+#' function for idframe objects
+#' 
+#' @return an idframe object
+#' 
+#' @examples
+#' data(cstr)
+#' cstrsub <- dataSlice(cstr,start=200,end=400) # extract between indices 200 and 400
+#' cstrTrain <- dataSlice(cstr,end=4500) # extract upto index 4500
+#' cstrTest <- dataSlice(cstr,start=6501) # extract from index 6501 till the end
+#' cstr_new <- dataSlice(cstr,freq=0.5) # resample data at half the original frequency  
+#' 
+#' @seealso \code{\link[stats]{window}}
+#' @export
+dataSlice <- function(data,start=NULL,end=NULL,freq=NULL){
+  # check if the class is correct
+  if(class(data)!='idframe')
+    stop("Not an idframe data")
+  
+  indexWindow <- function(y,start,end,freq){
+    Y <- matrix(y,ncol=ncol(y)); z <- as.vector(time(y))
+    Y <- window(Y,start=start,end=end,frequency=freq)
+    zw <- window(z,start=start,end=end,frequency=freq)
+    temp <- ts(Y,start=zw[1],end=tail(zw,n=1),deltat=diff(zw)[1])
+    colnames(temp) <- colnames(y)
+    temp
+  }
+  if(nOutputSeries(data)!=0)
+    outputData(data) <- indexWindow(outputData(data),start,end,freq)
+  
+  if(nInputSeries(data)!=0)
+    inputData(data) <- indexWindow(inputData(data),start,end,freq)
+  
+  return(data)
+}
\ No newline at end of file
diff --git a/R/rarx.R b/R/rarx.R
new file mode 100644
index 0000000..7f48888
--- /dev/null
+++ b/R/rarx.R
@@ -0,0 +1,79 @@
+#' Estimate parameters of ARX recursively
+#' 
+#' Estimates the parameters of a single-output ARX model of the 
+#' specified order from data using the recursive weighted least-squares
+#' algorithm.
+#' 
+#' @param x an object of class \code{idframe}
+#' @param order Specification of the orders: the three integer components 
+#' (na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) and 
+#' the input-output delay
+#' @param lambda Forgetting factor(Default=\code{0.95}) 
+#' 
+#' @return 
+#' A list containing the following objects
+#' \describe{
+#'  \item{theta}{Estimated parameters of the model. The \eqn{k^{th}} 
+#'  row contains the parameters associated with the \eqn{k^{th}} 
+#'  sample. Each row in \code{theta} has the following format: \cr
+#'  theta[i,:]=[a1,a2,...,ana,b1,...bnb]
+#'  }
+#'  \item{yhat}{Predicted value of the output, according to the 
+#'  current model - parameters based on all past data}
+#' }
+#' 
+#' @references
+#' Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+#' Theory and Practice}, CRC Press, Boca Raton. Section 25.1.3
+#' 
+#' Lennart Ljung (1999), \emph{System Identification: Theory for the User}, 
+#' 2nd Edition, Prentice Hall, New York. Section 11.2
+#' @examples 
+#' Gp1 <- idpoly(c(1,-0.9,0.2),2,ioDelay=2,noiseVar = 0.1)
+#' Gp2 <- idpoly(c(1,-1.2,0.35),2.5,ioDelay=2,noiseVar = 0.1)
+#' uk = idinput(2044,'prbs',c(0,1/4)); N = length(uk);
+#' N1 = round(0.35*N); N2 = round(0.4*N); N3 = N-N1-N2;
+#' yk1 <- sim(Gp1,uk[1:N1],addNoise = TRUE)
+#' yk2 <- sim(Gp2,uk[N1+1:N2],addNoise = TRUE)
+#' yk3 <- sim(Gp1,uk[N1+N2+1:N3],addNoise = TRUE)
+#' yk <- c(yk1,yk2,yk3)
+#' z <- idframe(yk,uk,1)
+#' g(theta,yhat) %=% rarx(z,c(2,1,2))
+#' 
+#' @export
+rarx <- function(x,order=c(1,1,1),lambda=0.95){
+  y <- outputData(x); u <- inputData(x)
+  N <- dim(y)[1]
+  na <- order[1];nb <- order[2]; nk <- order[3]
+  nb1 <- nb+nk-1 ; n <- max(na,nb1); df <- N-na-nb
+  
+  yout <- apply(y,2,padZeros,n=n)
+  uout <- apply(u,2,padZeros,n=n)
+  
+  uindex <- nk:nb1
+  if(na!=0) yindex <- 1:na
+  
+  reg <- function(i) {
+    # regressor
+    temp <- numeric(0)
+    if(na!=0) temp <- c(temp,-yout[i-yindex,])
+    phi <- c(temp,uout[i-uindex,])
+    phi
+  }
+  
+  # R0 <- reg(n+1)%*%t(reg(n+1))
+  # Plast <- solve(R0)
+  Plast <- 10^4*diag(na+nb)
+  theta <- matrix(0,N+1,na+nb)
+  yhat <- y
+  
+  for(i in 1:N){
+    temp <- reg(n+i)
+    yhat[i,] <- t(temp)%*%t(theta[i,,drop=FALSE])
+    eps_i <- y[i,,drop=FALSE] - yhat[i,,drop=FALSE]
+    kappa_i <- Plast%*%temp/(lambda+t(temp)%*%Plast%*%temp)[1]
+    theta[i+1,] <- t(t(theta[i,,drop=F])+eps_i[1]*kappa_i)
+    Plast <- (diag(na+nb)-kappa_i%*%t(temp))%*%Plast/lambda
+  }
+  list(theta=theta[1+1:N,],yhat=yhat)
+}
\ No newline at end of file
diff --git a/R/readData.R b/R/readData.R
new file mode 100644
index 0000000..33c1dcb
--- /dev/null
+++ b/R/readData.R
@@ -0,0 +1,73 @@
+#' Data input into a idframe object
+#' 
+#' Read the contents of a data.frame/matrix into a \code{idframe} object.
+#' 
+#' @param data a \code{data.frame} object
+#' @param ninputs the number of input columns. (Default: 0)
+#' @param Ts sampling interval (Default: 1)
+#' @param unit Time Unit (Default: "seconds")
+#' 
+#' 
+#' @return an idframe object
+#' @examples
+#' data(cstrData)
+#' data <- read.idframe(cstrData,ninputs=1,Ts= 1,unit="minutes")
+#' 
+#' @export 
+read.idframe <- function(data,ninputs=NULL,Ts = 1,
+                         unit=c("seconds","minutes","hours",
+                                "days")[1]){
+    outIndex <- 1:dim(data)[2]; inputs <- NULL
+    if(ninputs!=0){
+      inputs <- data[,1:ninputs,drop=F]
+      outIndex <- seq(ninputs+1,dim(data)[2],by=1)
+    }
+    outputs <- data[,outIndex,drop=F]
+    
+    out <- idframe(output=outputs,input=inputs,Ts=Ts,unit=unit)
+  
+  return(out)
+}
+
+#' Read the contents of a table-formatted file
+#' 
+#' Read the contents of an file in table format into a \code{idframe} object.
+#' 
+#' @param file the path to the file to read
+#' @param sep the field separator character. Values on each line of the file are 
+#' separated by this character. (Default: \code{","})
+#' @param header a logical value indicating whether the first row corresponding to 
+#' the first element of the rowIndex vector contains the names of the variables. 
+#' (Default: \code{TRUE})
+#' @param ninputs the number of input columns. (Default: 0)
+#' @param Ts sampling interval (Default: 1)
+#' @param unit Time Unit (Default: "seconds")
+#' @param ... additional arguments to be passed to the \code{\link[utils]{read.table}} function
+#' 
+#' @details
+#' 
+#' The \code{read.table.idframe} function uses the \code{\link[utils]{read.table}} function, 
+#' provided by the \pkg{utils} package, to read data from a table-formatted file and then calls the 
+#' \code{\link{read.idframe}} function to read the data into a idframe object
+#' 
+#' @return an idframe object
+#' @examples
+#' dataMatrix <- data.frame(matrix(rnorm(1000),ncol=5))
+#' colnames(dataMatrix) <- c("u1","u2","y1","y2","y3")
+#' write.csv(dataMatrix,file="test.csv",row.names=FALSE)
+#'  
+#' data <- read.table.idframe("test.csv",ninputs=2,unit="minutes")
+#' 
+#' @seealso  \code{\link[utils]{read.table}}
+#' @export 
+read.table.idframe <- function(file,header=TRUE,sep=",",ninputs=0,
+             Ts = 1,unit=c("seconds","minutes","hours",
+                           "days")[1],...){
+  
+  # Read from file (default: csv file)
+  dat <- read.table(file=file,header=header,sep=sep,...)
+  
+  # read from dataframe and return idframe object
+  out <- read.idframe(dat,ninputs=ninputs,Ts = Ts,unit=unit)
+  return(out)
+}
\ No newline at end of file
diff --git a/R/sim.R b/R/sim.R
new file mode 100644
index 0000000..2efa74e
--- /dev/null
+++ b/R/sim.R
@@ -0,0 +1,87 @@
+#' Simulate response of dynamic system
+#' 
+#' Simulate the response of a system to a given input
+#' 
+#' @param model the linear system to simulate
+#' @param input a vector/matrix containing the input
+#' @param addNoise logical variable indicating whether to add noise to the 
+#' response model. (Default: \code{FALSE})
+#' @param innov an optional times series of innovations. If not supplied (specified
+#' as \code{NULL}), gaussian white noise is generated, with the variance specified in 
+#' the model (Property: \code{noiseVar})
+#' @param seed integer indicating the seed value of the random number generator.
+#' Useful for reproducibility purposes.
+#' 
+#' @return
+#' a vector containing the simulated output
+#' 
+#' @details
+#' The routine is currently built only for SISO systems. Future versions will
+#' include support for MIMO systems.
+#' 
+#' @examples
+#' # ARX Model
+#' u <- idinput(300,"rgs")
+#' model <- idpoly(A=c(1,-1.5,0.7),B=c(0.8,-0.25),ioDelay=1,
+#' noiseVar=0.1)
+#' y <- sim(model,u,addNoise=TRUE)
+#' 
+#' @export
+sim <- function(model,input,addNoise = F,innov=NULL,seed=NULL) UseMethod("sim")
+
+#' @export
+sim.default <- function(model,input,addNoise = F,innov=NULL,seed=NULL){
+  print("The sim method is not developed for the current class of the object")
+}
+
+#' @import signal polynom
+#' @export
+sim.idpoly <- function(model,input,addNoise = F,innov=NULL,seed=NULL){
+  B <- c(rep(0,model$ioDelay),model$B)
+  Gden <- as.numeric(polynomial(model$A)*polynomial(model$F1))
+  G <- signal::Arma(b=B,a=Gden)
+  ufk <- signal::filter(G,input)
+  yk <- as.numeric(ufk)
+  
+  if(addNoise){
+    n <- ifelse(is.numeric(input),length(input),dim(input)[1])
+    if(!is.null(innov)){
+      ek <- innov
+    } else{
+      if(!is.null(seed)) set.seed(seed)
+      ek <- rnorm(n,sd=sqrt(model$noiseVar))
+    }
+    den1 <- as.numeric(polynomial(model$A)*polynomial(model$D))
+    filt1 <- Arma(b=model$C,a=den1)
+    vk <- signal::filter(filt1,ek)
+    
+    yk <- yk + as.numeric(vk)
+  }
+  
+  return(yk) 
+}
+
+# sim_arx <- function(model,input,ek){
+#   na <- length(model$A) - 1; nk <- model$ioDelay; 
+#   nb <- length(model$B) - 1; nb1 <- nb+nk
+#   n <- max(na,nb1)
+#   coef <- matrix(c(model$A[-1],model$B),nrow=na+(nb+1))
+#   
+#   if(class(input)=="idframe"){
+#     uk <- input$input[,1,drop=T]
+#   } else if(class(input) %in% c("matrix","data.frame")){
+#     uk <- input[,1,drop=T]
+#   } else if(is.numeric(input)){
+#     uk <- input
+#   }
+#   
+#   y <- rep(0,length(input)+n)
+#   u <- c(rep(0,n),uk)
+#   
+#   for(i in n+1:length(uk)){
+#     if(nk==0) v <- u[i-0:nb] else v <- u[i-nk:nb1]
+#     reg <- matrix(c(-(y[i-1:na]),v),ncol=na+(nb+1))
+#     y[i] <- reg%*%coef + ek[i]
+#   }
+#   return(y[n+1:length(uk)])
+# }
\ No newline at end of file
diff --git a/R/sysid.R b/R/sysid.R
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/R/sysid.R
diff --git a/R/util.R b/R/util.R
new file mode 100644
index 0000000..4667f9a
--- /dev/null
+++ b/R/util.R
@@ -0,0 +1,53 @@
+#' Multiple assignment operator
+#' 
+#' Assign multiple variables from a list or function return object
+#' 
+#' @param l the variables to be assigned
+#' @param r the list or function-return object
+#' 
+#' @aliases g
+#' 
+#' @export
+'%=%' = function(l,r) UseMethod('%=%')
+
+# Binary Operator
+#' @export
+'%=%.lbunch' = function(l,r) {
+  Envir = as.environment(-1)
+  
+#   if (length(r) > length(l))
+#     warning("RHS has more args than LHS. Only first", length(l), "used.")
+  
+  if (length(l) > length(r))  {
+    warning("LHS has more args than RHS. RHS will be repeated.")
+    r <- extendToMatch(r, l)
+  }
+  
+  for (II in 1:length(l)) {
+    do.call('<-', list(l[[II]], r[[II]]), envir=Envir)
+  }
+}
+
+# Used if LHS is larger than RHS
+extendToMatch <- function(source, destin) {
+  s <- length(source)
+  d <- length(destin)
+  
+  # Assume that destin is a length when it is a single number and source is not
+  if(d==1 && s>1 && !is.null(as.numeric(destin)))
+    d <- destin
+  
+  dif <- d - s
+  if (dif > 0) {
+    source <- rep(source, ceiling(d/s))[1:d]
+  }
+  return (source)
+}
+
+# Grouping the left hand side
+#' @export
+g = function(...) {
+  List = as.list(substitute(list(...)))[-1L]
+  class(List) = 'lbunch'
+  return(List)
+}
\ No newline at end of file
diff --git a/README.md b/README.md
index a852e29..46e5cd4 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,3 @@
-# SysID-R-code
-R code for system identification developed by the team at IIT Madras
+# sysid - System Identification in R
+
+The package provides functions for constructing mathematical models of dynamic systems from measured input-output data. It lets you create and use models of dynamic systems not easily modeled from first principles or specifications. You can use time-domain and frequency-domain input-output data to identify discrete-time input-output models, transfer functions and state space models.
diff --git a/data/armaxsim.rda b/data/armaxsim.rda
new file mode 100644
index 0000000..e7c4501
--- /dev/null
+++ b/data/armaxsim.rda
diff --git a/data/arxsim.rda b/data/arxsim.rda
new file mode 100644
index 0000000..0faed61
--- /dev/null
+++ b/data/arxsim.rda
diff --git a/data/bjsim.rda b/data/bjsim.rda
new file mode 100644
index 0000000..3c49c27
--- /dev/null
+++ b/data/bjsim.rda
diff --git a/data/cstr.rda b/data/cstr.rda
new file mode 100644
index 0000000..d518282
--- /dev/null
+++ b/data/cstr.rda
diff --git a/data/cstrData.rda b/data/cstrData.rda
new file mode 100644
index 0000000..01389f0
--- /dev/null
+++ b/data/cstrData.rda
diff --git a/data/cstr_mis.rda b/data/cstr_mis.rda
new file mode 100644
index 0000000..00851d2
--- /dev/null
+++ b/data/cstr_mis.rda
diff --git a/data/frd.rda b/data/frd.rda
new file mode 100644
index 0000000..7f847f3
--- /dev/null
+++ b/data/frd.rda
diff --git a/data/oesim.rda b/data/oesim.rda
new file mode 100644
index 0000000..f0e1c92
--- /dev/null
+++ b/data/oesim.rda
diff --git a/man/armax.Rd b/man/armax.Rd
new file mode 100644
index 0000000..18310d7
--- /dev/null
+++ b/man/armax.Rd
@@ -0,0 +1,71 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/estpoly.R
+\name{armax}
+\alias{armax}
+\title{Estimate ARMAX Models}
+\usage{
+armax(x, order = c(0, 1, 1, 0), init_sys = NULL, intNoise = FALSE,
+  options = optimOptions())
+}
+\arguments{
+\item{x}{an object of class \code{idframe}}
+
+\item{order}{Specification of the orders: the four integer components 
+(na,nb,nc,nk) are the order of polynolnomial A, order of polynomial B 
++ 1, order of the polynomial C,and the input-output delay respectively}
+
+\item{init_sys}{Linear polynomial model that configures the initial parameterization.
+Must be an ARMAX model. Overrules the \code{order} argument}
+
+\item{intNoise}{Logical variable indicating whether to add integrators in
+the noise channel (Default=\code{FALSE})}
+
+\item{options}{Estimation Options, setup using \code{\link{optimOptions}}}
+}
+\value{
+An object of class \code{estpoly} containing the following elements:
+ \item{sys}{an \code{idpoly} object containing the 
+   fitted ARMAX coefficients}
+ \item{fitted.values}{the predicted response}
+ \item{residuals}{the residuals}
+ \item{input}{the input data used}
+ \item{call}{the matched call}
+ \item{stats}{A list containing the following fields: \cr
+     \code{vcov} - the covariance matrix of the fitted coefficients \cr
+     \code{sigma} - the standard deviation of the innovations}
+ \item{options}{Option set used for estimation. If no 
+   custom options were configured, this is a set of default options}
+ \item{termination}{Termination conditions for the iterative
+    search used for prediction error minimization:
+     \code{WhyStop} - Reason for termination \cr
+     \code{iter} - Number of Iterations \cr
+     \code{iter} - Number of Function Evaluations }
+}
+\description{
+Fit an ARMAX model of the specified order given the input-output data
+}
+\details{
+SISO ARMAX models are of the form 
+\deqn{
+   y[k] + a_1 y[k-1] + \ldots + a_{na} y[k-na] = b_{nk} u[k-nk] + 
+   \ldots + b_{nk+nb} u[k-nk-nb] + c_{1} e[k-1] + \ldots c_{nc} e[k-nc]
+   + e[k] 
+}
+The function estimates the coefficients using non-linear least squares 
+(Levenberg-Marquardt Algorithm)
+\cr
+The data is expected to have no offsets or trends. They can be removed 
+using the \code{\link{detrend}} function.
+}
+\examples{
+data(armaxsim)
+z <- dataSlice(armaxsim,end=1533) # training set
+mod_armax <- armax(z,c(1,2,1,2))
+mod_armax
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+Theory and Practice}, CRC Press, Boca Raton. Sections 14.4.1, 21.6.2
+}
+
diff --git a/man/armaxsim.Rd b/man/armaxsim.Rd
new file mode 100644
index 0000000..69f41f1
--- /dev/null
+++ b/man/armaxsim.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{armaxsim}
+\alias{armaxsim}
+\title{Data simulated from an ARMAX model}
+\format{an \code{idframe} object with 2555 samples, one input and one
+output}
+\usage{
+armaxsim
+}
+\description{
+This dataset contains 2555 samples simulated from the following ARMAX model:
+\deqn{
+ y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + 
+ \frac{1-0.3q^{-1}}{1 - 0.5q^{-1}} e[k]
+}
+}
+\details{
+The model is simulated with a 2555 samples long full-band PRBS input.
+The noise variance is set to 0.1
+}
+\keyword{datasets}
+
diff --git a/man/arx.Rd b/man/arx.Rd
new file mode 100644
index 0000000..f55db17
--- /dev/null
+++ b/man/arx.Rd
@@ -0,0 +1,70 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/estpoly.R
+\name{arx}
+\alias{arx}
+\title{Estimate ARX Models}
+\usage{
+arx(x, order = c(1, 1, 1), lambda = 0.1, intNoise = FALSE, fixed = NULL)
+}
+\arguments{
+\item{x}{an object of class \code{idframe}}
+
+\item{order}{Specification of the orders: the three integer components 
+(na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) and 
+the input-output delay}
+
+\item{lambda}{Regularization parameter(Default=\code{0.1})}
+
+\item{intNoise}{Logical variable indicating whether to add integrators in
+the noise channel (Default=\code{FALSE})}
+
+\item{fixed}{list containing fixed parameters. If supplied, only \code{NA} entries 
+will be varied. Specified as a list of two vectors, each containing the parameters
+of polynomials A and B respectively.}
+}
+\value{
+An object of class \code{estpoly} containing the following elements:
+ \item{sys}{an \code{idpoly} object containing the 
+   fitted ARX coefficients}
+ \item{fitted.values}{the predicted response}
+ \item{residuals}{the residuals}
+ \item{input}{the input data used}
+ \item{call}{the matched call}
+ \item{stats}{A list containing the following fields: \cr
+     \code{vcov} - the covariance matrix of the fitted coefficients \cr
+     \code{sigma} - the standard deviation of the innovations\cr
+     \code{df} - the residual degrees of freedom}
+}
+\description{
+Fit an ARX model of the specified order given the input-output data
+}
+\details{
+SISO ARX models are of the form 
+\deqn{
+   y[k] + a_1 y[k-1] + \ldots + a_{na} y[k-na] = b_{nk} u[k-nk] + 
+   \ldots + b_{nk+nb} u[k-nk-nb] + e[k] 
+}
+The function estimates the coefficients using linear least squares (with
+regularization).
+\cr
+The data is expected to have no offsets or trends. They can be removed 
+using the \code{\link{detrend}} function. 
+\cr
+To estimate finite impulse response(\code{FIR}) models, specify the first
+order to be zero.
+}
+\examples{
+data(arxsim)
+mod_arx <- arx(arxsim,c(1,2,2))
+mod_arx
+plot(mod_arx) # plot the predicted and actual responses
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+Theory and Practice}, CRC Press, Boca Raton. Section 21.6.1
+
+Lennart Ljung (1999), \emph{System Identification: Theory for the User}, 
+2nd Edition, Prentice Hall, New York. Section 10.1
+}
+
diff --git a/man/arxsim.Rd b/man/arxsim.Rd
new file mode 100644
index 0000000..4972ff6
--- /dev/null
+++ b/man/arxsim.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{arxsim}
+\alias{arxsim}
+\title{Data simulated from an ARX model}
+\format{an \code{idframe} object with 2555 samples, one input and one
+output}
+\usage{
+arxsim
+}
+\description{
+This dataset contains 2555 samples simulated from the following ARX model:
+\deqn{
+ y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + 
+ \frac{1}{1 - 0.5q^{-1}} e[k]
+}
+}
+\details{
+The model is simulated with a 2555 samples long full-band PRBS input.
+The noise variance is set to 0.1
+}
+\keyword{datasets}
+
diff --git a/man/bj.Rd b/man/bj.Rd
new file mode 100644
index 0000000..b6a7cfc
--- /dev/null
+++ b/man/bj.Rd
@@ -0,0 +1,86 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/estpoly.R
+\name{bj}
+\alias{bj}
+\title{Estimate Box-Jenkins Models}
+\usage{
+bj(z, order = c(1, 1, 1, 1, 0), init_sys = NULL, options = optimOptions())
+}
+\arguments{
+\item{z}{an \code{idframe} object containing the data}
+
+\item{order}{Specification of the orders: the five integer components 
+(nb,nc,nd,nf,nk) are order of polynomial B + 1, order of the polynomial C,
+order of the polynomial D, order of the polynomial F, and the 
+input-output delay respectively}
+
+\item{init_sys}{Linear polynomial model that configures the initial parameterization.
+Must be a BJ model. Overrules the \code{order} argument}
+
+\item{options}{Estimation Options, setup using 
+\code{\link{optimOptions}}}
+}
+\value{
+An object of class \code{estpoly} containing the following elements:
+ \item{sys}{an \code{idpoly} object containing the 
+   fitted BJ coefficients}
+ \item{fitted.values}{the predicted response}
+ \item{residuals}{the residuals}
+ \item{input}{the input data used}
+ \item{call}{the matched call}
+ \item{stats}{A list containing the following fields: \cr
+     \code{vcov} - the covariance matrix of the fitted coefficients \cr
+     \code{sigma} - the standard deviation of the innovations}
+ \item{options}{Option set used for estimation. If no 
+   custom options were configured, this is a set of default options}
+ \item{termination}{Termination conditions for the iterative
+    search used for prediction error minimization:
+     \code{WhyStop} - Reason for termination \cr
+     \code{iter} - Number of Iterations \cr
+     \code{iter} - Number of Function Evaluations }
+}
+\description{
+Fit a box-jenkins model of the specified order from input-output data
+}
+\details{
+SISO BJ models are of the form 
+\deqn{
+   y[k] = \frac{B(q^{-1})}{F(q^{-1})}u[k-nk] + 
+   \frac{C(q^{-1})}{D(q^{-1})} e[k]
+}
+The orders of Box-Jenkins model are defined as follows:
+\deqn{
+   B(q^{-1}) = b_1 + b_2q^{-1} + \ldots + b_{nb} q^{-nb+1}
+}
+
+\deqn{
+   C(q^{-1}) = 1 + c_1q^{-1} + \ldots + c_{nc} q^{-nc}
+}
+
+\deqn{
+   D(q^{-1}) = 1 + d_1q^{-1} + \ldots + d_{nd} q^{-nd}
+}
+\deqn{
+   F(q^{-1}) = 1 + f_1q^{-1} + \ldots + f_{nf} q^{-nf}
+}
+
+The function estimates the coefficients using non-linear least squares 
+(Levenberg-Marquardt Algorithm)
+\cr
+The data is expected to have no offsets or trends. They can be removed 
+using the \code{\link{detrend}} function.
+}
+\examples{
+data(bjsim)
+z <- dataSlice(bjsim,end=1500) # training set
+mod_bj <- bj(z,c(2,1,1,1,2))
+mod_bj 
+residplot(mod_bj) # residual plots
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+Theory and Practice}, CRC Press, Boca Raton. Sections 14.4.1, 17.5.2, 
+21.6.3
+}
+
diff --git a/man/bjsim.Rd b/man/bjsim.Rd
new file mode 100644
index 0000000..d599860
--- /dev/null
+++ b/man/bjsim.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{bjsim}
+\alias{bjsim}
+\title{Data simulated from an BJ model}
+\format{an \code{idframe} object with 2046 samples, one input and one
+output}
+\usage{
+bjsim
+}
+\description{
+This dataset contains 2046 samples simulated from the following BJ model:
+\deqn{
+ y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + 
+ \frac{1+0.2q^{-1}}{1 - 0.3q^{-1}} e[k]
+}
+}
+\details{
+The model is simulated with a 2046 samples long full-band PRBS input.
+The noise variance is set to 0.1
+}
+\keyword{datasets}
+
diff --git a/man/compare.Rd b/man/compare.Rd
new file mode 100644
index 0000000..3700f0a
--- /dev/null
+++ b/man/compare.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/predict.R
+\name{compare}
+\alias{compare}
+\title{Compare the measured output and the predicted output(s)}
+\usage{
+compare(data, nahead = 1, ...)
+}
+\arguments{
+\item{data}{validation data in the form of an \code{idframe} object}
+
+\item{nahead}{number of steps ahead at which to predict (Default:1). For infinite-
+step ahead predictions, supply \code{Inf}.}
+
+\item{\ldots}{models whose predictions are to be compared}
+}
+\description{
+Plots the output predictions of model(s) superimposed over validation data, 
+data, for comparison.
+}
+\examples{
+data(arxsim)
+mod1 <- arx(arxsim,c(1,2,2))
+mod2 <- oe(arxsim,c(2,1,1))
+compare(arxsim,nahead=Inf,mod1,mod2)
+
+}
+\seealso{
+\code{\link{predict.estpoly}} for obtaining model predictions
+}
+
diff --git a/man/cstr.Rd b/man/cstr.Rd
new file mode 100644
index 0000000..66f31f5
--- /dev/null
+++ b/man/cstr.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{cstr}
+\alias{cstr}
+\title{Continuous stirred tank reactor data (idframe)}
+\format{an \code{idframe} object with 7500 samples, one input and two
+outputs}
+\usage{
+cstr
+}
+\description{
+The Process is a model of a Continuous Stirring Tank Reactor, 
+where the reaction is exothermic and the concentration is 
+controlled by regulating the coolant flow.
+\cr
+}
+\details{
+Inputs: q, Coolant Flow l/min
+Outputs:
+\describe{
+\item{Ca}{Concentration mol/l}
+\item{T}{Temperature Kelvin}}
+}
+\keyword{datasets}
+
diff --git a/man/cstrData.Rd b/man/cstrData.Rd
new file mode 100644
index 0000000..7fc7be9
--- /dev/null
+++ b/man/cstrData.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{cstrData}
+\alias{cstrData}
+\title{Continuous stirred tank reactor data (data.frame)}
+\format{an \code{data.frame} object with 7500 rows and three columns:
+q, Ca and T}
+\source{
+\url{ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/cstr.dat.gz}
+}
+\usage{
+cstrData
+}
+\description{
+The Process is a model of a Continuous Stirring Tank Reactor, 
+where the reaction is exothermic and the concentration is 
+controlled by regulating the coolant flow.
+\cr
+}
+\details{
+Inputs: q, Coolant Flow l/min
+Outputs:
+\describe{
+\item{Ca}{Concentration mol/l}
+\item{T}{Temperature Kelvin}}
+}
+\keyword{datasets}
+
diff --git a/man/cstr_mis.Rd b/man/cstr_mis.Rd
new file mode 100644
index 0000000..15f9433
--- /dev/null
+++ b/man/cstr_mis.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{cstr_mis}
+\alias{cstr_mis}
+\title{Continuous stirred tank reactor data with missing values}
+\format{an \code{idframe} object with 7500 samples, one input and two
+outputs}
+\usage{
+cstr_mis
+}
+\description{
+This dataset is derived from the \code{cstr} dataset with few samples
+containing missing values, in one or all variables. It is used to 
+demonstrate the capabilities of the \code{misdata} routine.
+}
+\seealso{
+\code{\link{cstr}}, \code{\link{misdata}}
+}
+\keyword{datasets}
+
diff --git a/man/dataSlice.Rd b/man/dataSlice.Rd
new file mode 100644
index 0000000..f0c1cd6
--- /dev/null
+++ b/man/dataSlice.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/preprocess.R
+\name{dataSlice}
+\alias{dataSlice}
+\title{Subset or Resample idframe data}
+\usage{
+dataSlice(data, start = NULL, end = NULL, freq = NULL)
+}
+\arguments{
+\item{data}{an object of class \code{idframe}}
+
+\item{start}{the start index}
+
+\item{end}{the end index}
+
+\item{freq}{fraction of the original frequency at which the series
+to be sampled.}
+}
+\value{
+an idframe object
+}
+\description{
+\code{dataSlice} is a subsetting method for objects of class \code{idframe}. It 
+extracts the subset of the object \code{data} observed between indices \code{start}
+and \code{end}.  If a frequency is specified, the series is then re-sampled at the 
+new frequency.
+}
+\details{
+The dataSlice function extends the \code{\link[stats]{window}} 
+function for idframe objects
+}
+\examples{
+data(cstr)
+cstrsub <- dataSlice(cstr,start=200,end=400) # extract between indices 200 and 400
+cstrTrain <- dataSlice(cstr,end=4500) # extract upto index 4500
+cstrTest <- dataSlice(cstr,start=6501) # extract from index 6501 till the end
+cstr_new <- dataSlice(cstr,freq=0.5) # resample data at half the original frequency  
+
+}
+\seealso{
+\code{\link[stats]{window}}
+}
+
diff --git a/man/detrend.Rd b/man/detrend.Rd
new file mode 100644
index 0000000..2303b5b
--- /dev/null
+++ b/man/detrend.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/preprocess.R
+\name{detrend}
+\alias{detrend}
+\alias{trInfo}
+\title{Remove offsets and linear trends}
+\usage{
+detrend(x, type = 0)
+}
+\arguments{
+\item{x}{an object of class \code{idframe}}
+
+\item{type}{argument indicating the type of trend to be removed (Default=\code{0})
+\itemize{
+  \item type=\code{0}: Subtracts mean value from each signal
+  \item type=\code{1}: Subtracts a linear trend (least-squres fit)
+  \item type=\code{trInfo} object: Subtracts a trend specified by the object
+}}
+}
+\value{
+A list containing two objects: the detrended data and the trend information
+}
+\description{
+Removes offsets or trends from data
+}
+\details{
+\code{R} by default doesn't allow return of multiple objects. The \code{\%=\%}
+operator and \code{g} function in this package facillitate this behaviour. See 
+the examples section for more information.
+}
+\examples{
+data(cstr)
+datatrain <- dataSlice(cstr,end=4500)
+datatest <- dataSlice(cstr,4501)
+g(Ztrain,tr) \%=\% detrend(datatrain) # Remove means
+g(Ztest) \%=\% detrend(datatest,tr)
+
+}
+\seealso{
+\code{\link[stats]{lm}}
+}
+
diff --git a/man/estpoly.Rd b/man/estpoly.Rd
new file mode 100644
index 0000000..13a5556
--- /dev/null
+++ b/man/estpoly.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/estpoly.R
+\name{estpoly}
+\alias{estpoly}
+\title{Estimated polynomial object}
+\usage{
+estpoly(sys, fitted.values, residuals, options = NULL, call, stats,
+  termination = NULL, input)
+}
+\arguments{
+\item{sys}{an \code{idpoly} object containing the estimated polynomial
+coefficients}
+
+\item{fitted.values}{1-step ahead predictions on the training dataset}
+
+\item{residuals}{1-step ahead prediction errors}
+
+\item{options}{optimization specification ser used (applicable for non-linear least
+squares)}
+
+\item{call}{the matched call}
+
+\item{stats}{a list containing estimation statistics}
+
+\item{termination}{termination criteria for optimization}
+
+\item{input}{input signal of the training data-set}
+}
+\description{
+Estimated discrete-time polynomial model returned from an estimation
+routine.
+}
+\details{
+Do not use \code{estpoly} for directly specifing an input-output polynomial model.
+\code{\link{idpoly}} is to be used instead
+}
+
diff --git a/man/etfe.Rd b/man/etfe.Rd
new file mode 100644
index 0000000..662ff4a
--- /dev/null
+++ b/man/etfe.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nonparam.R
+\name{etfe}
+\alias{etfe}
+\title{Estimate empirical transfer function}
+\usage{
+etfe(data, n = 128)
+}
+\arguments{
+\item{data}{an object of class \code{idframe}}
+
+\item{n}{frequency spacing (Default: \code{128})}
+}
+\value{
+an \code{idfrd} object containing the estimated frequency response
+}
+\description{
+Estimates the emperical transfer function from the data by taking the 
+ratio of the fourier transforms of the output and the input variables
+}
+\examples{
+data(arxsim)
+frf <- etfe(arxsim)
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+Theory and Practice}, CRC Press, Boca Raton. Sections 5.3 and 20.4.2
+}
+\seealso{
+\code{\link[stats]{fft}}
+}
+
diff --git a/man/fitch.Rd b/man/fitch.Rd
new file mode 100644
index 0000000..1715797
--- /dev/null
+++ b/man/fitch.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/estpoly.R
+\name{fitch}
+\alias{fitch}
+\title{Fit Characteristics}
+\usage{
+fitch(x)
+}
+\arguments{
+\item{x}{the estimated model}
+}
+\value{
+A list containing the following elements
+
+\item{MSE}{Mean Square Error measure of how well the response of the model fits
+the estimation data}
+\item{FPE}{Final Prediction Error}
+\item{FitPer}{Normalized root mean squared error (NRMSE) measure of how well the 
+response of the model fits the estimation data, expressed as a percentage.}
+\item{AIC}{Raw Akaike Information Citeria (AIC) measure of model quality}
+\item{AICc}{Small sample-size corrected AIC}
+\item{nAIC}{Normalized AIC}
+\item{BIC}{Bayesian Information Criteria (BIC)}
+}
+\description{
+Returns quantitative assessment of the estimated model as a list
+}
+
diff --git a/man/frd.Rd b/man/frd.Rd
new file mode 100644
index 0000000..c122fed
--- /dev/null
+++ b/man/frd.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{frd}
+\alias{frd}
+\title{Frequency response data}
+\format{an \code{idfrd} object with response at 128 frequency points}
+\usage{
+frd
+}
+\description{
+This dataset contains frequency response data of an unknown SISO system.
+}
+\keyword{datasets}
+
diff --git a/man/getcov.Rd b/man/getcov.Rd
new file mode 100644
index 0000000..8ee2537
--- /dev/null
+++ b/man/getcov.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/estUtil.R
+\name{getcov}
+\alias{getcov}
+\title{Parameter covariance of the identified model}
+\usage{
+getcov(sys)
+}
+\arguments{
+\item{sys}{a linear, identified parametric model}
+}
+\description{
+Obtain the parameter covariance matrix of the linear, identified 
+parametric model
+}
+
diff --git a/man/grapes-equals-grapes.Rd b/man/grapes-equals-grapes.Rd
new file mode 100644
index 0000000..fbc3231
--- /dev/null
+++ b/man/grapes-equals-grapes.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/util.R
+\name{\%=\%}
+\alias{\%=\%}
+\alias{g}
+\title{Multiple assignment operator}
+\usage{
+l \%=\% r
+}
+\arguments{
+\item{l}{the variables to be assigned}
+
+\item{r}{the list or function-return object}
+}
+\description{
+Assign multiple variables from a list or function return object
+}
+
diff --git a/man/idframe.Rd b/man/idframe.Rd
new file mode 100644
index 0000000..ce00ba0
--- /dev/null
+++ b/man/idframe.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/idframe.R
+\name{idframe}
+\alias{idframe}
+\title{S3 class for storing input-output data.}
+\usage{
+idframe(output, input = NULL, Ts = 1, start = 0, end = NULL,
+  unit = c("seconds", "minutes", "hours", "days")[1])
+}
+\arguments{
+\item{output}{dataframe/matrix/vector containing the outputs}
+
+\item{input}{dataframe/matrix/vector containing the inputs}
+
+\item{Ts}{sampling interval (Default: 1)}
+
+\item{start}{Time of the first observation}
+
+\item{end}{Time of the last observation Optional Argument}
+
+\item{unit}{Time unit (Default: "seconds")}
+}
+\value{
+an idframe object
+}
+\description{
+\code{idframe} is an S3 class for storing and manipulating input-ouput data. It supports discrete time and frequency domain data.
+}
+\examples{
+
+dataMatrix <- matrix(rnorm(1000),ncol=5) 
+data <- idframe(output=dataMatrix[,3:5],input=dataMatrix[,1:2],Ts=1)
+
+}
+\seealso{
+\code{\link{plot.idframe}}, the plot method for idframe objects
+}
+
diff --git a/man/idfrd.Rd b/man/idfrd.Rd
new file mode 100644
index 0000000..1d3f96e
--- /dev/null
+++ b/man/idfrd.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/idframe.R
+\name{idfrd}
+\alias{idfrd}
+\title{S3 class constructor for storing frequency response data}
+\usage{
+idfrd(respData, freq, Ts, spec = NULL, covData = NULL, noiseCov = NULL)
+}
+\arguments{
+\item{respData}{frequency response data. For SISO systems, supply a 
+vector of frequency response values. For MIMO systems with Ny 
+outputs and Nu inputs, supply an array of size c(Ny,Nu,Nw).}
+
+\item{freq}{frequency points of the response}
+
+\item{Ts}{sampling time of data}
+
+\item{spec}{power spectra and cross spectra of the system 
+output disturbances (noise). Supply an array of size (Ny,Ny,Nw)}
+
+\item{covData}{response data covariance matrices. Supply an array
+of size (Ny,Nu,Nw,2,2). covData[ky,ku,kw,,] is the covariance matrix
+of respData[ky,ku,kw]}
+
+\item{noiseCov}{power spectra variance. Supply an array of 
+size (Ny,Ny,Nw)}
+}
+\value{
+an idfrd object
+}
+\description{
+S3 class constructor for storing frequency response data
+}
+\seealso{
+\code{\link{plot.idfrd}} for generating bode plots, 
+\code{\link{spa}} and \code{\link{etfe}} for estimating the 
+frequency response given input/output data
+}
+
diff --git a/man/idinput.Rd b/man/idinput.Rd
new file mode 100644
index 0000000..01f0b7c
--- /dev/null
+++ b/man/idinput.Rd
@@ -0,0 +1,36 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/idinput.R
+\name{idinput}
+\alias{idinput}
+\title{function to generate input singals (rgs/rbs/prbs/sine)}
+\usage{
+idinput(n, type = "rgs", band = c(0, 1), levels = c(-1, 1))
+}
+\arguments{
+\item{n}{integer length of the input singal to be generated}
+
+\item{type}{the type of input signal to be generated. 
+'rgs' - generates random gaussian signal
+'rbs' - generates random binary signal
+'prbs' - generates pseudorandom binary signal
+'sine' - generates a signal that is a sum of sinusoids
+
+Default value is type='rgs'}
+
+\item{band}{determines the frequency content of the signal. 
+For type='rbs'/'sine'/,  band = [wlow,whigh]
+which specifies the lower and the upper bound of the passband frequencies(expressed as fractions of Nyquist frequency). Default is c(0,1)
+For type='prbs', band=[0,B]
+where B is such that the singal is constant over 1/B (clock period). Default is c(0,1)}
+
+\item{levels}{row vector defining the input level. It is of the form 
+levels=c(minu, maxu)
+For 'rbs','prbs', 'sine', the generated signal always between minu and maxu.
+For 'rgs', minu=mean value of signal minus one standard deviation and maxu=mean value of signal plus one standard deviation
+
+Default value is levels=c(-1,1)}
+}
+\description{
+\code{idinput} is a function for generating input signals (rgs/rbs/prbs/sine) for identification purposes
+}
+
diff --git a/man/idpoly.Rd b/man/idpoly.Rd
new file mode 100644
index 0000000..f95eca2
--- /dev/null
+++ b/man/idpoly.Rd
@@ -0,0 +1,54 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/poly.R
+\name{idpoly}
+\alias{idpoly}
+\title{Polynomial model with identifiable parameters}
+\usage{
+idpoly(A = 1, B = 1, C = 1, D = 1, F1 = 1, ioDelay = 0, Ts = 1,
+  noiseVar = 1, intNoise = F, unit = c("seconds", "minutes", "hours",
+  "days")[1])
+}
+\arguments{
+\item{A}{autoregressive coefficients}
+
+\item{B, F1}{coefficients of the numerator and denominator respectively
+of the deterministic model between the input and output}
+
+\item{C, D}{coefficients of the numerator and denominator respectively
+of the stochastic model}
+
+\item{ioDelay}{the delay in the input-output channel}
+
+\item{Ts}{sampling interval}
+
+\item{noiseVar}{variance of the white noise source (Default=\code{1})}
+
+\item{intNoise}{Logical variable indicating presence or absence of integrator
+in the noise channel (Default=\code{FALSE})}
+
+\item{unit}{time unit (Default=\code{"seconds"})}
+}
+\description{
+Creates a polynomial model with identifiable coefficients
+}
+\details{
+Discrete-time polynomials are of the form
+\deqn{
+ A(q^{-1}) y[k] = \frac{B(q^{-1})}{F1(q^{-1})} u[k] + 
+ \frac{C(q^{-1})}{D(q^{-1})} e[k] 
+}
+}
+\examples{
+# define output-error model
+mod_oe <- idpoly(B=c(0.6,-0.2),F1=c(1,-0.5),ioDelay = 2,Ts=0.1,
+noiseVar = 0.1)
+
+# define box-jenkins model with unit variance
+B <- c(0.6,-0.2)
+C <- c(1,-0.3)
+D <- c(1,1.5,0.7)
+F1 <- c(1,-0.5)
+mod_bj <- idpoly(1,B,C,D,F1,ioDelay=1)
+
+}
+
diff --git a/man/impulseest.Rd b/man/impulseest.Rd
new file mode 100644
index 0000000..4247a0f
--- /dev/null
+++ b/man/impulseest.Rd
@@ -0,0 +1,46 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nonparam.R
+\name{impulseest}
+\alias{impulseest}
+\title{Estimate Impulse Response Coefficients}
+\usage{
+impulseest(x, M = 30, K = NULL, regul = F, lambda = 1)
+}
+\arguments{
+\item{x}{an object of class \code{idframe}}
+
+\item{M}{Order of the FIR Model (Default:\code{30})}
+
+\item{K}{Transport delay in the estimated impulse response 
+(Default:NULL)}
+
+\item{regul}{Parameter indicating whether regularization should be 
+used. (Default:\code{FALSE})}
+
+\item{lambda}{The value of the regularization parameter. Valid only if
+\code{regul=TRUE}. (Default:\code{1})}
+}
+\description{
+\code{impulseest} is used to estimate impulse response coefficients from 
+the data
+}
+\details{
+The IR Coefficients are estimated using linear least squares. Future 
+Versions will provide support for multivariate data.
+}
+\examples{
+uk <- rnorm(1000,1)
+yk <- filter (uk,c(0.9,-0.4),method="recursive") + rnorm(1000,1)
+data <- idframe(output=data.frame(yk),input=data.frame(uk))
+fit <- impulseest(data)
+impulseplot(fit)
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+Theory and Practice}, CRC Press, Boca Raton. Sections 17.4.11 and 20.2
+}
+\seealso{
+\code{\link{step}}
+}
+
diff --git a/man/impulseplot.Rd b/man/impulseplot.Rd
new file mode 100644
index 0000000..0e7ced3
--- /dev/null
+++ b/man/impulseplot.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nonparam.R
+\name{impulseplot}
+\alias{impulseplot}
+\title{Impulse Response Plots}
+\usage{
+impulseplot(model, sd = 2)
+}
+\arguments{
+\item{model}{an object of class \code{impulseest}}
+
+\item{sd}{standard deviation of the confidence region (Default: \code{2})}
+}
+\description{
+Plots the estimated IR coefficients along with the significance limits
+at each lag.
+}
+\seealso{
+\code{\link{impulseest}},\code{\link{step}}
+}
+
diff --git a/man/inputData.Rd b/man/inputData.Rd
new file mode 100644
index 0000000..5164aa8
--- /dev/null
+++ b/man/inputData.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ioNamesData.R
+\name{inputData}
+\alias{inputData}
+\alias{inputData.idframe}
+\alias{outputData}
+\alias{outputData.idframe}
+\title{Output or Input-data}
+\usage{
+inputData(x, series)
+}
+\arguments{
+\item{x}{\code{idframe} object}
+
+\item{series}{the indices to extract}
+}
+\description{
+Extract output-data or input-data in idframe objects
+}
+
diff --git a/man/inputNames-set.Rd b/man/inputNames-set.Rd
new file mode 100644
index 0000000..7584a55
--- /dev/null
+++ b/man/inputNames-set.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ioNamesData.R
+\name{inputNames<-}
+\alias{inputNames}
+\alias{inputNames<-}
+\alias{inputNames<-.idframe}
+\alias{outputNames}
+\alias{outputNames<-}
+\alias{outputNames<-.idframe}
+\title{Extract or set series' names}
+\usage{
+inputNames(x) <- value
+}
+\arguments{
+\item{x}{\code{idframe} object}
+
+\item{value}{vector of strings}
+}
+\description{
+Extract or set names of series in input or output
+}
+
diff --git a/man/iv.Rd b/man/iv.Rd
new file mode 100644
index 0000000..91a4340
--- /dev/null
+++ b/man/iv.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/iv.R
+\name{iv}
+\alias{iv}
+\title{ARX model estimation using instrumental variable method}
+\usage{
+iv(z, order = c(0, 1, 0), x = NULL)
+}
+\arguments{
+\item{z}{an idframe object containing the data}
+
+\item{order}{Specification of the orders: the three integer components 
+(na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) 
+and the input-output delay}
+
+\item{x}{instrument variable matrix. x must be of the same size as the output 
+data. (Default: \code{NULL})}
+}
+\value{
+An object of class \code{estpoly} containing the following elements:
+ \item{sys}{an \code{idpoly} object containing the 
+   fitted ARX coefficients}
+ \item{fitted.values}{the predicted response}
+ \item{residuals}{the residuals}
+ \item{input}{the input data used}
+ \item{call}{the matched call}
+ \item{stats}{A list containing the following fields: \cr
+     \code{vcov} - the covariance matrix of the fitted coefficients \cr
+     \code{sigma} - the standard deviation of the innovations\cr
+     \code{df} - the residual degrees of freedom}
+}
+\description{
+Estimates an ARX model of the specified order from input-output data using
+the instrument variable method. If arbitrary instruments are not supplied 
+by the user, the instruments are generated using the arx routine
+}
+\examples{
+data(arxsim)
+mod_iv <- iv(arxsim,c(2,1,1))
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+Theory and Practice}, CRC Press, Boca Raton. Sections 21.7.1, 21.7.2
+
+Lennart Ljung (1999), \emph{System Identification: Theory for the User}, 
+2nd Edition, Prentice Hall, New York. Section 7.6
+}
+\seealso{
+\code{\link{arx}}, \code{\link{iv4}}
+}
+
diff --git a/man/iv4.Rd b/man/iv4.Rd
new file mode 100644
index 0000000..77b794c
--- /dev/null
+++ b/man/iv4.Rd
@@ -0,0 +1,55 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/iv.R
+\name{iv4}
+\alias{iv4}
+\title{ARX model estimation using four-stage instrumental variable method}
+\usage{
+iv4(z, order = c(0, 1, 0))
+}
+\arguments{
+\item{z}{an idframe object containing the data}
+
+\item{order}{Specification of the orders: the three integer components 
+(na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) 
+and the input-output delay}
+}
+\value{
+An object of class \code{estpoly} containing the following elements:
+ \item{sys}{an \code{idpoly} object containing the 
+   fitted ARX coefficients}
+ \item{fitted.values}{the predicted response}
+ \item{residuals}{the residuals}
+ \item{input}{the input data used}
+ \item{call}{the matched call}
+ \item{stats}{A list containing the following fields: \cr
+     \code{vcov} - the covariance matrix of the fitted coefficients \cr
+     \code{sigma} - the standard deviation of the innovations\cr
+     \code{df} - the residual degrees of freedom}
+}
+\description{
+Estimates an ARX model of the specified order from input-output data using
+the instrument variable method. The estimation algorithm is insensitive to 
+the color of the noise term.
+}
+\details{
+Estimation is performed in 4 stages. The first stage uses the arx function. The resulting model generates the 
+instruments for a second-stage IV estimate. The residuals obtained from this model are modeled using a sufficently
+high-order AR model. At the fourth stage, the input-output data is filtered through this AR model and then subjected 
+to the IV function with the same instrument filters as in the second stage.
+}
+\examples{
+mod_dgp <- idpoly(A=c(1,-0.5),B=c(0.6,-.2),C=c(1,0.6),ioDelay = 2,noiseVar = 0.1)
+u <- idinput(400,"prbs")
+y <- sim(mod_dgp,u,addNoise=TRUE)
+z <- idframe(y,u)
+mod_iv4 <- iv4(z,c(1,2,2))
+
+}
+\references{
+Lennart Ljung (1999), \emph{System Identification: Theory for the User}, 
+2nd Edition, Prentice Hall, New York. Section 15.3
+}
+\seealso{
+\code{\link{arx}}, \code{\link{iv4}}
+}
+
diff --git a/man/misdata.Rd b/man/misdata.Rd
new file mode 100644
index 0000000..a6d4df8
--- /dev/null
+++ b/man/misdata.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/preprocess.R
+\name{misdata}
+\alias{misdata}
+\title{Replace Missing Data by Interpolation}
+\usage{
+misdata(data)
+}
+\arguments{
+\item{data}{an object of class \code{idframe}}
+}
+\value{
+data (an idframe object) with missing data replaced.
+}
+\description{
+Function for replacing missing values with interpolated ones. This is an
+extension of the \code{na.approx} function from the \code{zoo} package.
+The missing data is indicated using the value \emph{NA}.
+}
+\examples{
+data(cstr_mis)
+summary(cstr_mis) # finding out the number of NAs
+cstr <- misdata(cstr_mis)
+
+}
+\seealso{
+\code{\link[zoo]{na.approx}}
+}
+
diff --git a/man/nInputSeries.Rd b/man/nInputSeries.Rd
new file mode 100644
index 0000000..7f7feed
--- /dev/null
+++ b/man/nInputSeries.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ioNamesData.R
+\name{nInputSeries}
+\alias{nInputSeries}
+\alias{nOutputSeries}
+\title{Number of series in input or output}
+\usage{
+nInputSeries(data)
+}
+\arguments{
+\item{data}{\code{idframe} object}
+}
+\description{
+Number of series in input or output in a idframe object
+}
+
diff --git a/man/oe.Rd b/man/oe.Rd
new file mode 100644
index 0000000..f5a56bc
--- /dev/null
+++ b/man/oe.Rd
@@ -0,0 +1,70 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/estpoly.R
+\name{oe}
+\alias{oe}
+\title{Estimate Output-Error Models}
+\usage{
+oe(x, order = c(1, 1, 0), init_sys = NULL, options = optimOptions())
+}
+\arguments{
+\item{x}{an object of class \code{idframe}}
+
+\item{order}{Specification of the orders: the four integer components 
+(nb,nf,nk) are order of polynomial B + 1, order of the polynomial F,
+and the input-output delay respectively}
+
+\item{init_sys}{Linear polynomial model that configures the initial parameterization.
+Must be an OE model. Overrules the \code{order} argument}
+
+\item{options}{Estimation Options, setup using 
+\code{\link{optimOptions}}}
+}
+\value{
+An object of class \code{estpoly} containing the following elements:
+ \item{sys}{an \code{idpoly} object containing the 
+   fitted OE coefficients}
+ \item{fitted.values}{the predicted response}
+ \item{residuals}{the residuals}
+ \item{input}{the input data used}
+ \item{call}{the matched call}
+ \item{stats}{A list containing the following fields: \cr
+     \code{vcov} - the covariance matrix of the fitted coefficients \cr
+     \code{sigma} - the standard deviation of the innovations}
+ \item{options}{Option set used for estimation. If no 
+   custom options were configured, this is a set of default options}
+ \item{termination}{Termination conditions for the iterative
+    search used for prediction error minimization:
+     \code{WhyStop} - Reason for termination \cr
+     \code{iter} - Number of Iterations \cr
+     \code{iter} - Number of Function Evaluations }
+}
+\description{
+Fit an output-error model of the specified order given the input-output data
+}
+\details{
+SISO OE models are of the form 
+\deqn{
+   y[k] + f_1 y[k-1] + \ldots + f_{nf} y[k-nf] = b_{nk} u[k-nk] + 
+   \ldots + b_{nk+nb} u[k-nk-nb] + f_{1} e[k-1] + \ldots f_{nf} e[k-nf]
+   + e[k] 
+}
+The function estimates the coefficients using non-linear least squares 
+(Levenberg-Marquardt Algorithm)
+\cr
+The data is expected to have no offsets or trends. They can be removed 
+using the \code{\link{detrend}} function.
+}
+\examples{
+data(oesim)
+z <- dataSlice(oesim,end=1533) # training set
+mod_oe <- oe(z,c(2,1,2))
+mod_oe 
+plot(mod_oe) # plot the predicted and actual responses
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+Theory and Practice}, CRC Press, Boca Raton. Sections 14.4.1, 17.5.2, 
+21.6.3
+}
+
diff --git a/man/oesim.Rd b/man/oesim.Rd
new file mode 100644
index 0000000..7ecf0a7
--- /dev/null
+++ b/man/oesim.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{oesim}
+\alias{oesim}
+\title{Data simulated from an OE model}
+\format{an \code{idframe} object with 2555 samples, one input and one
+output}
+\usage{
+oesim
+}
+\description{
+This dataset contains 2555 samples simulated from the following OE model:
+\deqn{
+ y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + e[k]
+}
+}
+\details{
+The model is simulated with a 2555 samples long full-band PRBS input.
+The noise variance is set to 0.1
+}
+\keyword{datasets}
+
diff --git a/man/optimOptions.Rd b/man/optimOptions.Rd
new file mode 100644
index 0000000..1d77f64
--- /dev/null
+++ b/man/optimOptions.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/estUtil.R
+\name{optimOptions}
+\alias{optimOptions}
+\title{Create optimization options}
+\usage{
+optimOptions(tol = 0.01, maxIter = 20, LMinit = 0.01, LMstep = 2,
+  display = c("off", "on")[1])
+}
+\arguments{
+\item{tol}{Minimum 2-norm of the gradient (Default: \code{1e-2})}
+
+\item{maxIter}{Maximum number of iterations to be performed}
+
+\item{LMinit}{Starting value of search-direction length 
+in the Levenberg-Marquardt method (Default: \code{0.01})}
+
+\item{LMstep}{Size of the Levenberg-Marquardt step (Default: \code{2})}
+
+\item{display}{Argument whether to display iteration details or not
+(Default: \code{"off"})}
+}
+\description{
+Specify optimization options that are to be passed to the 
+numerical estimation routines
+}
+
diff --git a/man/plot.idframe.Rd b/man/plot.idframe.Rd
new file mode 100644
index 0000000..2e7948d
--- /dev/null
+++ b/man/plot.idframe.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/idframe.R
+\name{plot.idframe}
+\alias{plot.idframe}
+\title{Plotting idframe objects}
+\usage{
+\method{plot}{idframe}(x, col = "steelblue", lwd = 1, main = NULL,
+  size = 12, ...)
+}
+\arguments{
+\item{x}{an \code{idframe} object}
+
+\item{col}{line color, to be passed to plot.(Default=\code{"steelblue"})}
+
+\item{lwd}{line width, in millimeters(Default=\code{1})}
+
+\item{main}{the plot title. (Default = \code{NULL})}
+
+\item{size}{text size (Default = \code{12})}
+
+\item{\ldots}{additional arguments}
+}
+\description{
+Plotting method for objects inherting from class \code{idframe}
+}
+\examples{
+data(cstr)
+plot(cstr,col="blue")
+
+}
+
diff --git a/man/plot.idfrd.Rd b/man/plot.idfrd.Rd
new file mode 100644
index 0000000..435f6d7
--- /dev/null
+++ b/man/plot.idfrd.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/idframe.R
+\name{plot.idfrd}
+\alias{plot.idfrd}
+\title{Plotting idfrd objects}
+\usage{
+\method{plot}{idfrd}(x, col = "steelblue", lwd = 1, ...)
+}
+\arguments{
+\item{x}{An object of class \code{idframe}}
+
+\item{col}{a specification for the line colour (Default : \code{"
+steelblue"})}
+
+\item{lwd}{the line width, a positive number, defaulting to 1}
+
+\item{\ldots}{additional arguments}
+}
+\description{
+Generates the bode plot of the given frequency response data. It uses the
+ggplot2 plotting engine
+}
+\examples{
+data(frd)
+plot(frd)
+
+}
+\seealso{
+\code{\link[ggplot2]{ggplot}}
+}
+
diff --git a/man/predict.estpoly.Rd b/man/predict.estpoly.Rd
new file mode 100644
index 0000000..53371e6
--- /dev/null
+++ b/man/predict.estpoly.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/predict.R
+\name{predict.estpoly}
+\alias{predict.estpoly}
+\title{Predictions of identified model}
+\usage{
+\method{predict}{estpoly}(object, newdata = NULL, nahead = 1, ...)
+}
+\arguments{
+\item{object}{\code{estpoly} object containing the identified model}
+
+\item{newdata}{optional dataset to be used for predictions. If not supplied, 
+predictions are made on the training set.}
+
+\item{nahead}{number of steps ahead at which to predict (Default:1). For infinite-
+step ahead predictions or pure simulation, supply \code{Inf}.}
+
+\item{\ldots}{other arguments}
+}
+\value{
+Time-series containing the predictions
+}
+\description{
+Predicts the output of an identified model (\code{estpoly}) object K steps ahead.
+}
+\examples{
+data(arxsim)
+mod1 <- oe(arxsim,c(2,1,1))
+Yhat <- predict(mod1,arxsim) #  1-step ahead predictions 
+Yhat_2 <- predict(mod1,arxsim,nahead=2) # 2-step ahead predictions
+Yhat_inf <- predict(mod1,arxsim,nahead=Inf) # Infinite-step ahead predictions
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: Theory 
+and Practice}, CRC Press, Boca Raton. Chapter 18
+}
+
diff --git a/man/rarx.Rd b/man/rarx.Rd
new file mode 100644
index 0000000..1d0bac4
--- /dev/null
+++ b/man/rarx.Rd
@@ -0,0 +1,55 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rarx.R
+\name{rarx}
+\alias{rarx}
+\title{Estimate parameters of ARX recursively}
+\usage{
+rarx(x, order = c(1, 1, 1), lambda = 0.95)
+}
+\arguments{
+\item{x}{an object of class \code{idframe}}
+
+\item{order}{Specification of the orders: the three integer components 
+(na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) and 
+the input-output delay}
+
+\item{lambda}{Forgetting factor(Default=\code{0.95})}
+}
+\value{
+A list containing the following objects
+\describe{
+ \item{theta}{Estimated parameters of the model. The \eqn{k^{th}} 
+ row contains the parameters associated with the \eqn{k^{th}} 
+ sample. Each row in \code{theta} has the following format: \cr
+ theta[i,:]=[a1,a2,...,ana,b1,...bnb]
+ }
+ \item{yhat}{Predicted value of the output, according to the 
+ current model - parameters based on all past data}
+}
+}
+\description{
+Estimates the parameters of a single-output ARX model of the 
+specified order from data using the recursive weighted least-squares
+algorithm.
+}
+\examples{
+Gp1 <- idpoly(c(1,-0.9,0.2),2,ioDelay=2,noiseVar = 0.1)
+Gp2 <- idpoly(c(1,-1.2,0.35),2.5,ioDelay=2,noiseVar = 0.1)
+uk = idinput(2044,'prbs',c(0,1/4)); N = length(uk);
+N1 = round(0.35*N); N2 = round(0.4*N); N3 = N-N1-N2;
+yk1 <- sim(Gp1,uk[1:N1],addNoise = TRUE)
+yk2 <- sim(Gp2,uk[N1+1:N2],addNoise = TRUE)
+yk3 <- sim(Gp1,uk[N1+N2+1:N3],addNoise = TRUE)
+yk <- c(yk1,yk2,yk3)
+z <- idframe(yk,uk,1)
+g(theta,yhat) \%=\% rarx(z,c(2,1,2))
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+Theory and Practice}, CRC Press, Boca Raton. Section 25.1.3
+
+Lennart Ljung (1999), \emph{System Identification: Theory for the User}, 
+2nd Edition, Prentice Hall, New York. Section 11.2
+}
+
diff --git a/man/read.idframe.Rd b/man/read.idframe.Rd
new file mode 100644
index 0000000..9d821b8
--- /dev/null
+++ b/man/read.idframe.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/readData.R
+\name{read.idframe}
+\alias{read.idframe}
+\title{Data input into a idframe object}
+\usage{
+read.idframe(data, ninputs = NULL, Ts = 1, unit = c("seconds", "minutes",
+  "hours", "days")[1])
+}
+\arguments{
+\item{data}{a \code{data.frame} object}
+
+\item{ninputs}{the number of input columns. (Default: 0)}
+
+\item{Ts}{sampling interval (Default: 1)}
+
+\item{unit}{Time Unit (Default: "seconds")}
+}
+\value{
+an idframe object
+}
+\description{
+Read the contents of a data.frame/matrix into a \code{idframe} object.
+}
+\examples{
+data(cstrData)
+data <- read.idframe(cstrData,ninputs=1,Ts= 1,unit="minutes")
+
+}
+
diff --git a/man/read.table.idframe.Rd b/man/read.table.idframe.Rd
new file mode 100644
index 0000000..ca80a4c
--- /dev/null
+++ b/man/read.table.idframe.Rd
@@ -0,0 +1,50 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/readData.R
+\name{read.table.idframe}
+\alias{read.table.idframe}
+\title{Read the contents of a table-formatted file}
+\usage{
+read.table.idframe(file, header = TRUE, sep = ",", ninputs = 0, Ts = 1,
+  unit = c("seconds", "minutes", "hours", "days")[1], ...)
+}
+\arguments{
+\item{file}{the path to the file to read}
+
+\item{header}{a logical value indicating whether the first row corresponding to 
+the first element of the rowIndex vector contains the names of the variables. 
+(Default: \code{TRUE})}
+
+\item{sep}{the field separator character. Values on each line of the file are 
+separated by this character. (Default: \code{","})}
+
+\item{ninputs}{the number of input columns. (Default: 0)}
+
+\item{Ts}{sampling interval (Default: 1)}
+
+\item{unit}{Time Unit (Default: "seconds")}
+
+\item{...}{additional arguments to be passed to the \code{\link[utils]{read.table}} function}
+}
+\value{
+an idframe object
+}
+\description{
+Read the contents of an file in table format into a \code{idframe} object.
+}
+\details{
+The \code{read.table.idframe} function uses the \code{\link[utils]{read.table}} function, 
+provided by the \pkg{utils} package, to read data from a table-formatted file and then calls the 
+\code{\link{read.idframe}} function to read the data into a idframe object
+}
+\examples{
+dataMatrix <- data.frame(matrix(rnorm(1000),ncol=5))
+colnames(dataMatrix) <- c("u1","u2","y1","y2","y3")
+write.csv(dataMatrix,file="test.csv",row.names=FALSE)
+ 
+data <- read.table.idframe("test.csv",ninputs=2,unit="minutes")
+
+}
+\seealso{
+\code{\link[utils]{read.table}}
+}
+
diff --git a/man/residplot.Rd b/man/residplot.Rd
new file mode 100644
index 0000000..b9509c9
--- /dev/null
+++ b/man/residplot.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/estpoly.R
+\name{residplot}
+\alias{residplot}
+\title{Plot residual characteristics}
+\usage{
+residplot(model, newdata = NULL)
+}
+\arguments{
+\item{model}{estimated polynomial model}
+
+\item{newdata}{an optional dataset on which predictions are to be computed. If
+not supplied, predictions are computed on the training dataset.}
+}
+\description{
+Computes the 1-step ahead prediction errors (residuals) for an estimated polynomial
+model, and plots auto-correlation of the residuals and the 
+cross-correlation of the residuals with the input signals.
+}
+
diff --git a/man/sim.Rd b/man/sim.Rd
new file mode 100644
index 0000000..a40ba29
--- /dev/null
+++ b/man/sim.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sim.R
+\name{sim}
+\alias{sim}
+\title{Simulate response of dynamic system}
+\usage{
+sim(model, input, addNoise = F, innov = NULL, seed = NULL)
+}
+\arguments{
+\item{model}{the linear system to simulate}
+
+\item{input}{a vector/matrix containing the input}
+
+\item{addNoise}{logical variable indicating whether to add noise to the 
+response model. (Default: \code{FALSE})}
+
+\item{innov}{an optional times series of innovations. If not supplied (specified
+as \code{NULL}), gaussian white noise is generated, with the variance specified in 
+the model (Property: \code{noiseVar})}
+
+\item{seed}{integer indicating the seed value of the random number generator.
+Useful for reproducibility purposes.}
+}
+\value{
+a vector containing the simulated output
+}
+\description{
+Simulate the response of a system to a given input
+}
+\details{
+The routine is currently built only for SISO systems. Future versions will
+include support for MIMO systems.
+}
+\examples{
+# ARX Model
+u <- idinput(300,"rgs")
+model <- idpoly(A=c(1,-1.5,0.7),B=c(0.8,-0.25),ioDelay=1,
+noiseVar=0.1)
+y <- sim(model,u,addNoise=TRUE)
+
+}
+
diff --git a/man/spa.Rd b/man/spa.Rd
new file mode 100644
index 0000000..306189f
--- /dev/null
+++ b/man/spa.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nonparam.R
+\name{spa}
+\alias{spa}
+\title{Estimate frequency response}
+\usage{
+spa(x, winsize = NULL, freq = NULL)
+}
+\arguments{
+\item{x}{an \code{idframe} object}
+
+\item{winsize}{lag size of the Hanning window (Default: \code{min
+(length(x)/10,30)})}
+
+\item{freq}{frequency points at which the response is evaluated
+(Default: \code{seq(1,128)/128*pi/Ts})}
+}
+\value{
+an \code{idfrd} object containing the estimated frequency response 
+and the noise spectrum
+}
+\description{
+Estimates frequency response and noise spectrum from data with 
+fixed resolution using spectral analysis
+}
+\examples{
+data(arxsim)
+frf <- spa(arxsim)
+
+}
+\references{
+Arun K. Tangirala (2015), \emph{Principles of System Identification: 
+Theory and Practice}, CRC Press, Boca Raton. Sections 16.5 and 20.4
+}
+
diff --git a/man/step.Rd b/man/step.Rd
new file mode 100644
index 0000000..a4c547d
--- /dev/null
+++ b/man/step.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nonparam.R
+\name{step}
+\alias{step}
+\title{Step Response Plots}
+\usage{
+step(model)
+}
+\arguments{
+\item{model}{an object of class \code{impulseest}}
+}
+\description{
+Plots the step response of a system, given the IR model
+}
+\examples{
+uk <- rnorm(1000,1)
+yk <- filter (uk,c(0.9,-0.4),method="recursive") + rnorm(1000,1)
+data <- idframe(output=data.frame(yk),input=data.frame(uk))
+fit <- impulseest(data)
+step(fit) 
+
+}
+\seealso{
+\code{\link{impulseest}}
+}
+
diff --git a/man/time.Rd b/man/time.Rd
new file mode 100644
index 0000000..91ba6c6
--- /dev/null
+++ b/man/time.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/idframe.R
+\name{time}
+\alias{deltat}
+\alias{frequency}
+\alias{time}
+\title{Sampling times of IO data
+ 
+\code{time} creates the vector of times at which data was sampled. \code{frequency} returns the number of damples per unit time and \code{deltat} the time-interval 
+between observations}
+\usage{
+time(x)
+}
+\arguments{
+\item{x}{a idframe object, or a univariate or multivariate time-series, or a vector or matrix}
+}
+\description{
+Sampling times of IO data
+ 
+\code{time} creates the vector of times at which data was sampled. \code{frequency} returns the number of damples per unit time and \code{deltat} the time-interval 
+between observations
+}
+