window.var_
function(y, X, fit, m, tvec, step = 1)
{
#
# PURPOSE:  compute window resampling variance estimates for
#           binary spatial time series
#
#
##
###  Let "m" be length of the window
##
#
#  m <- trunc( C * N^{d/(d+2)} )
#
#
##
###  Let N be the number of time points
##
#
    N <- length(unique(tvec))
    nobs <- length(y)   #
##
### compute number of sub-windows to consider
##
#
    n.windows <- trunc((N - m + 1)/step)    #
##
### preliminary calculations
##
#
    antilogit <- function(x)
    {
        exp(x)/(1 + exp(x))
    }
#
    p.fitted <- fitted(fit)
    v.fitted <- p.fitted * (1 - p.fitted)
    R <- (y - p.fitted)
    U.N <- t(X) %*% R
    H.N <- t(X) %*% ((v.fitted) * X)    #
##
### compute U U^T and the one-step U U^T 
##
# 
    p <- ncol(X)
    eUU <- matrix(0, p, p)
    eUU.1 <- matrix(0, p, p)
    mU <- matrix(0, p, 1)   #
    utimes <- unique(tvec)
    utimes <- utimes[order(utimes)]
    for(j in 1:n.windows) {
        start <- utimes[(j - 1) * step + 1]
        stop <- start + m - 1   #
        keep <- (tvec >= start) & (tvec <= stop)    #
        X.j <- X[keep,  ]
        U.j <- t(X.j) %*% R[keep]   #
        H.j <- t(X.j) %*% ((v.fitted[keep]) * X.j)
        beta.1 <- fit$coef + as.vector(solve(H.N - H.j) %*% (U.N - U.j))    #
        lp.j <- as.vector(X.j %*% beta.1)
        U.j.1 <- t(X.j) %*% (y[keep] - antilogit(lp.j)) #
        eUU <- eUU + ((N/m) * (U.j %*% t(U.j)))/n.windows
        eUU.1 <- eUU.1 + ((N/m) * (U.j.1 %*% t(U.j.1)))/n.windows
        mU <- mU + U.j.1/n.windows
    }
#
    I.N <- solve(H.N)   #
    center <- T
    if(center)
        eUU.1 <- eUU.1 - (N/m) * mU %*% t(mU)
    var.est.0 <- I.N %*% eUU %*% I.N
    var.est.1 <- I.N %*% eUU.1 %*% I.N  #
    out <- list(title = "Window resampling for estimating functions", 
        m = m, step = step, n.windows = 
        n.windows, independence = I.N, 
        window.0 = var.est.0, window.1 = var.est.1) #
    out
}