1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
#' @export
plot.estPoly <- function(model,newdata=NULL){
require(ggplot2)
require(reshape2)
if(is.null(newdata)){
ypred <- fitted(model)
yact <- fitted(model) + resid(model)
time <- model$time
} else{
if(class(newdata)!="idframe") stop("Only idframe objects allowed")
ypred <- sim(coef(model),newdata$input[,1,drop=F])
yact <- newdata$output[,1,drop=F]
time <- seq(from=newdata$t.start,to=newdata$t.end,by=newdata$Ts)
}
df <- data.frame(Predicted=ypred,Actual=yact,Time=time)
meltdf <- melt(df,id="Time")
ggplot(meltdf, aes(x = Time,y=value,colour=variable,group=variable)) +
geom_line() + theme_bw()
}
#' @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 - predict(model,newdata)
u <- newdata$input
}
acorr <- acf(e,plot = F); ccorr <- ccf(as.numeric(u),e,plot = F)
par(mfrow=c(2,1),mar=c(3,4,4,2))
plot(acorr,main="ACF of residuals")
plot(ccorr,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 data 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 and
#' the input-output delay
#'
#' @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
#' no regularization). Future versions may include regularization
#' parameters as well
#' \\
#' The data is expected to have no offsets or trends. They can be removed
#' using the \code{\link{detrend}} function.
#'
#' @return
#' An object with classes \code{estARX} and \code{estPoly}, containing
#' the following elements:
#'
#' \tabular{ll}{
#' \code{coefficients} \tab an \code{idpoly} object containing the
#' fitted coefficients \cr
#' \code{vcov} \tab the covariance matrix of the fitted coefficients\cr
#' \code{sigma} \tab the standard deviation of the innovations\cr
#' \code{df} \tab the residual degrees of freedom \cr
#' \code{fitted.values} \tab the predicted response \cr
#' \code{residuals} \tab the residuals \cr
#' \code{call} \tab the matched call \cr
#' \code{time} \tab the time of the data used \cr
#' \code{input} \tab the input data used
#' }
#'
#'
#' @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)
#' model <- estARX(data,c(2,1,1))
#' summary(model) # obtain estimates and their covariances
#' plot(model) # plot the predicted and actual responses
#'
#' @export
estARX <- function(data,order=c(0,1,0)){
y <- as.matrix(data$output)
u <- as.matrix(data$input); N <- dim(y)[1]
na <- order[1];nb <- order[2]; nk <- order[3]
nb1 <- nb+nk ; n <- max(na,nb1); df <- N - na - nb -nk
padZeros <- function(x,n) c(rep(0,n),x,rep(0,n))
yout <- apply(y,2,padZeros,n=n);
uout <- apply(u,2,padZeros,n=n);
reg <- function(i) cbind(-yout[i-1:na,],uout[i-nk:nb1])
X <- t(sapply(n+1:(N+n),reg))
Y <- yout[n+1:(N+n),,drop=F]
qx <- qr(X); coef <- qr.solve(qx,Y)
sigma2 <- sum((Y-X%*%coef)^2)/df
vcov <- sigma2 * chol2inv(qx$qr)
model <- arx(A = c(1,coef[1:na]),B = coef[na+1:nb1],ioDelay = nk)
time <- seq(from=data$t.start,to=data$t.end,by=data$Ts)
est <- list(coefficients = model,vcov = vcov, sigma = sqrt(sigma2),
df = df,fitted.values=(X%*%coef)[1:N,],
residuals=(Y-X%*%coef)[1:N,],call=match.call(),
time=time,input=u)
class(est) <- c("estARX","estPoly")
est
}
#' @export
predict.estARX <- function(model,newdata=NULL){
if(is.null(newdata)){
return(fitted(model))
} else{
return(sim(coef(model),as.numeric(newdata$input)))
}
}
#' @export
summary.estARX <- function(object)
{
coefs <- c(coef(object)$A[-1],coef(object)$B)
se <- sqrt(diag(object$vcov))
tval <- coefs / se
TAB <- cbind(Estimate = coefs,
StdErr = se,
t.value = tval,
p.value = 2*pt(-abs(tval), df=object$df))
na <- length(coef(object)$A) - 1; nk <- coef(object)$ioDelay;
nb <- length(coef(object)$B) - nk
rownames(TAB) <- rep("a",nrow(TAB))
for(i in 1:na) rownames(TAB)[i] <- paste("a",i,sep="")
for(j in (na+1):nrow(TAB)) {
rownames(TAB)[j] <- paste("b",j-na-1+nk,sep="")
}
res <- list(call=object$call,coefficients=TAB,sigma=object$sigma,
df=object$df)
class(res) <- "summary.estARX"
res
}
#' @export
print.summary.estARX <- function(object){
cat("Discrete-time ARX model: A(q^{-1})y[k] = B(q^{-1})u[k] + e[k] \n")
cat("Call: ");print(object$call);cat("\n\n")
print(coef(object))
cat(paste("\nsigma:",format(object$sigma,digits=4)))
cat(paste("\nDoF:",object$df))
}
|
