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
|
#' @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]
nb1 <- nb+nk-1 ; n <- max(na,nb1); df <- N-na-nb
if(is.null(x)){
# Initial Guess using ARX
mod_arx <- arx(z,order)
x <- matrix(sim(mod_arx$sys,u,sigma=0))
}
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))
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_iv1 <- iv(z,order)
A <- signal::Ma(c(1,mod_iv1$sys$A))
B <- signal::Ma(c(rep(0,nk),mod_iv1$sys$B))
# Step 3 (AR Modeling)
w <- matrix(as.numeric(signal::filter(A,y)) -
as.numeric(signal::filter(B,u)))
mod_ar <- ar(w,aic = F,order=na+nb)
L <- signal::Ma(c(1,-mod_ar$ar))
# Step 4
}
|
