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
|
#' 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)
}
|
