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
|
#' 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})
#'
#' @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.
#'
#' @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(z,c(2,1,1))
#'
#' @seealso arx
#'
#' @export
iv <- function(z,order=c(0,1,0),x=NULL){
y <- outputData(z); u <- inputData(z); N <- dim(y)[1]
na <- order[1];nb <- order[2]; nk <- order[3]
if(is.null(x)){
# Initial Guess using ARX
mod_arx <- arx(z,order)
x <- matrix(sim(mod_arx$sys,u))
}
ivcompute(y,u,x,na,nb,nk,n,N,z$unit)
}
ivcompute <- function(y,u,x,na,nb,nk,n,N,unit){
nb1 <- nb+nk-1 ; n <- max(na,nb1); df <- N-na-nb
padZeros <- function(x,n) c(rep(0,n),x,rep(0,n))
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))
# Estimator
lhs <- t(psi)%*%phi; lhsinv <- solve(lhs)
theta <- lhsinv%*%t(psi)%*%Y
# 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)
}
#' @export
iv4 <- function(z,order=c(0,1,0)){
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
# Steps 1-2
mod_arx <- arx(z,order)
x <- matrix(sim(mod_arx$sys,u))
mod_iv <- ivcompute(y,u,x,na,nb,nk,n,N,z$unit)
# Step 3
w <- resid(mod_iv)
mod_ar <- ar(w,aic = F,order=na+nb)
Lhat <- mod_ar$ar
# Step 4
# G2 <- signal::Arma(as.numeric(B),as.numeric(A))
# x2 <- predict(mod_iv1)
Lf <- function(x,L,...) matrix(as.numeric(stats::filter(x,L,...)))
filtered <- lapply(list(y,u,x),Lf,L=Lhat,method="convolution",
sides=1,circular = T)
yf <- filtered[[1]]; uf<- filtered[[2]]; xf <- filtered[[3]]
ivcompute(yf,uf,xf,na,nb,nk,n,N,z$unit)
}
|
