1 files changed, 22 insertions, 44 deletions

diff --git a/R/sim.R b/R/sim.R
index 834313a..8811210 100644
--- a/R/sim.R
+++ b/R/sim.R
@@ -2,8 +2,10 @@
 #' 
 #' Simulate the response of a system given the input
 #' 
-#' @param model the system model to simulate
+#' @param model the linear system to simulate
 #' @param input a vector/matrix containing the input
+#' @param innov an optional times series of innovations. If not provided,
+#' innovations are generated using the \code{rnorm} function
 #' @param sigma standard deviation of the innovations (Default= \code{0})
 #' @param seed integer indicating the seed value of the random number generator
 #' 
@@ -12,57 +14,39 @@
 #' 
 #' @details
 #' The routine is currently built only for SISO systems. Future Versions will
-#' include support for MIMO systems. Current support 
+#' include support for MIMO systems.
 #' 
-#' @seealso 
-#' \code{\link{sim.idpoly}} for simulating polynomial models
+#' @examples
+#' # ARX Model
+#' u <- rnorm(200,sd=1)
+#' model <- idpoly(A=c(1,-1.5,0.7),B=c(0.8,-0.25),ioDelay=1)
+#' y <- sim(model,u,sigma=0.1)
 #' 
 #' @export
-sim <- function(model,input,sigma=0,seed=NULL) UseMethod("sim")
+sim <- function(model,input,innov=NULL,sigma=0,seed=NULL) UseMethod("sim")
 
 #' @export
 sim.default <- function(model,input,sigma=0,seed=NULL){
   print("The sim method is not developed for the current class of the object")
 }
 
-#' Simulate from a Polynomial Model
-#' 
-#' Simulate the response of a system governed by a polynomial model
-#' , given the input
-#' 
-#' @param model an object of class \code{idpoly} containing the coefficients
-#' @param input a vector/matrix containing the input
-#' @param sigma standard deviation of the innovations (Default= \code{0})
-#' @param seed integer indicating the seed value of the random number generator
-#' 
-#' @return
-#' a vector containing the output
-#' 
-#' @details
-#' The routine is currently built only for SISO systems. Future Versions will
-#' include support for MIMO systems
-#' 
-#' @seealso 
-#' \code{\link{idpoly}} for defining polynomial models
-#' 
-#' @examples
-#' # ARX Model
-#' u <- rnorm(200,sd=1)
-#' model <- idpoly(A=c(1,-1.5,0.7),B=c(0.8,-0.25),ioDelay=1)
-#' y <- sim(model,u,sigma=0.1)
-#' 
 #' @import polynom
-#' 
 #' @export
-sim.idpoly <- function(model,input,sigma=0,seed=NULL){
+sim.idpoly <- function(model,input,innov=NULL,sigma=0,seed=NULL){
+  if(!is.null(innov)){
+    ek <- innov
+  } else{
+    if(!is.null(seed)) set.seed(seed)
+    ek <- rnorm(n,sd=sigma)
+  }
+  
   if(model$type=="arx"){
-    sim_arx(model,input,sigma,seed)
+    sim_arx(model,input,ek)
   } else{
     require(signal);require(polynom)
     
     n <- length(input)[1]
-    if(!is.null(seed)) set.seed(seed)
-    ek <- rnorm(n,sd=sigma)
+   
     den1 <- as.numeric(polynomial(model$A)*polynomial(model$D))
     filt1 <- Arma(b=model$C,a=den1)
     vk <- signal::filter(filt1,ek)
@@ -78,7 +62,7 @@ sim.idpoly <- function(model,input,sigma=0,seed=NULL){
   }
 }
 
-sim_arx <- function(model,input,sigma=0,seed=NULL){
+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)
@@ -95,16 +79,10 @@ sim_arx <- function(model,input,sigma=0,seed=NULL){
   y <- rep(0,length(input)+n)
   u <- c(rep(0,n),uk)
   
-  if(!is.null(seed)) set.seed(seed)
-  
-  ek <- rnorm(length(uk),sd=sigma)
-  # padLeftZeros <- function(x) c(rep(0,n),x)
-  # u <- apply(input,2,padLeftZeros)
-  
   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-n]
+    y[i] <- reg%*%coef + ek[i]
   }
   return(y[n+1:length(uk)])
 }
\ No newline at end of file