Nothing
R/gs2slshet.R
In sphet: Estimation of Spatial Autoregressive Models with and without Heteroskedastic Innovations
Defines functions
impacts.gstsls
gstslshet
Documented in
gstslshet
impacts.gstsls
#' GM estimation of a Cliff-Ord type model with Heteroskedastic Innovations
#' @aliases gstslshet
#' @name gstslshet
#'
#' @description Multi step GM/IV estimation of a linear Cliff and Ord -type of model
#' of the form:
#'
#' \deqn{y=\lambda W y + X \beta + u }
#' \deqn{u=\rho W u + e} with \deqn{e ~ N(0,\sigma^2_i) }
#'
#' The model allows for spatial lag in the dependent variable
#' and disturbances. The innovations in the disturbance process are assumed
#' heteroskedastic of an unknown form.
#'
#' @usage gstslshet(formula, data = list(), listw, na.action = na.fail,
#' zero.policy = NULL, initial.value = 0.2, abs.tol = 1e-20,
#' rel.tol = 1e-10, eps = 1e-5, inverse = T, sarar = T)
#'
#' @param formula a description of the model to be fit
#' @param data an object of class \link{data.frame}. An optional data frame containing the variables
#' in the model.
#' @param listw an object of class \code{listw} created for example by \code{nb2listw}
#' @param na.action a function which indicates what should happen when the data contains missing values.
#' See \link{lm} for details.
#' @param zero.policy See \code{lagsarlm} for details
#' @param initial.value The initial value for \eqn{\rho}. It can be either numeric (default is 0.2) or
#' set to \code{'SAR'}, in which case the optimization will start from the estimated coefficient of a regression of the 2SLS
#' residuals over their spatial lag (i.e. a spatial AR model)
#' @param abs.tol Absolute tolerance. See \link{nlminb} for details.
#' @param rel.tol Relative tolerance. See \link{nlminb} for details.
#' @param eps Tolerance level for the approximation. See Details.
#' @param inverse \code{TRUE}. If \code{FALSE}, an appoximated inverse is calculated. See Details.
#' @param sarar \code{TRUE}. If \code{FALSE}, a spatial error model is estimated.
#'
#' @details The procedure consists of two steps alternating GM and IV estimators. Each step consists of sub-steps.
#' In step one \eqn{\delta = [\beta',\lambda]'} is estimated by 2SLS. The 2SLS residuals are first employed
#' to obtain an initial (consistent but not efficient) GM estimator of \eqn{\rho} and then a consistent and efficient
#' estimator (involving the variance-covariance matrix of the limiting distribution of the normalized sample moments).
#' In step two, the spatial Cochrane-Orcutt transformed model is estimated by 2SLS. This corresponds to a GS2SLS procedure.
#' The GS2SLS residuals are used to obtain a consistent and efficient GM estimator for \eqn{\rho}.
#'
#' The initial value for the optimization in step 1b is taken to be \code{initial.value}. The initial value in step 1c is the
#' optimal parameter of step 1b. Finally, the initial value for the optimization of step 2b is the optimal parameter of step 1c.
#'
#' Internally, the object of class \code{listw} is transformed into a \link{Matrix}
#' using the function \link{listw2dgCMatrix}.
#'
#'
#' The expression of the estimated variance covariance matrix of the limiting
#' distribution of the normalized sample moments based on 2SLS residuals
#' involves the inversion of \eqn{I-\rho W'}.
#' When \code{inverse} is \code{FALSE}, the inverse is calculated using the approximation
#' \eqn{I +\rho W' + \rho^2 W'^2 + ...+ \rho^n W'^n}.
#' The powers considered depend on a condition.
#' The function will keep adding terms until the absolute value of the \code{sum} of all elements
#' of the matrix \eqn{\rho^i W^i} is greater than a fixed \eqn{\epsilon} (\code{eps}). By default \code{eps}
#' is set to 1e-5.
#'
#'
#' @return A list object of class \code{sphet}
#' \item{coefficients}{Generalized Spatial two stage least squares coefficient estimates of \eqn{\delta} and GM estimator for \eqn{\rho}. }
#' \item{var}{variance-covariance matrix of the estimated coefficients}
#' \item{s2}{GS2SLS residuals variance}
#' \item{residuals}{GS2SLS residuals}
#' \item{yhat}{difference between GS2SLS residuals and response variable}
#' \item{call}{the call used to create this object}
#' \item{model}{the model matrix of data}
#' \item{method}{\code{'gs2slshac'}}
#' \item{W}{Wald test for both \eqn{\rho} and \eqn{\lambda} are zero}
#'
#' @references
#' Arraiz, I. and Drukker, M.D. and Kelejian, H.H. and Prucha, I.R. (2007)
#' A spatial Cliff-Ord-type Model with Heteroskedastic Innovations: Small and Large Sample Results,
#' \emph{Department of Economics, University of Maryland}'
#'
#' Kelejian, H.H. and Prucha, I.R. (2007)
#'Specification and Estimation of Spatial Autoregressive Models with Autoregressive and Heteroskedastic Disturbances,
#' \emph{Journal of Econometrics}, forthcoming.
#'
#' Kelejian, H.H. and Prucha, I.R. (1999)
#' A Generalized Moments Estimator for the Autoregressive Parameter in a Spatial Model,
#' \emph{International Economic Review}, \bold{40}, pages 509--533.
#'
#' Kelejian, H.H. and Prucha, I.R. (1998)
#'A Generalized Spatial Two Stage Least Square Procedure for Estimating a Spatial Autoregressive
#' Model with Autoregressive Disturbances,
#' \emph{Journal of Real Estate Finance and Economics}, \bold{17}, pages 99--121.
#'
#' @seealso \code{\link{stslshac}}
#'
#' @author Gianfranco Piras \email{gpiras@mac.com}
#'
#' @examples
#' data(columbus, package = "spdep")
#' listw <- spdep::nb2listw(col.gal.nb)
#' res <- gstslshet(CRIME ~ HOVAL + INC, data = columbus, listw = listw)
#' summary(res)
#'
#' @keywords models
#' @export
gstslshet<-function(formula, data=list(),listw, na.action=na.fail, zero.policy=NULL, initial.value=0.2, abs.tol=1e-20, rel.tol=1e-10, eps=1e-5, inverse=T, sarar=T){
##functions that need to be sourced
#source("twostagels.R")
#source("utilities.R")
#source("listw2dgCMatrix.R")
#source("Omega.R")
if (is.null(zero.policy))
zero.policy <- get.ZeroPolicyOption()
stopifnot(is.logical(zero.policy))
#extract model objects
mt<-terms(formula,data=data)
mf<-lm(formula, data, na.action=na.action, method="model.frame")
na.act<-attr(mf,'na.action')
# record call
cl<- match.call()
#if(!inverse) warning("Approximated inverse could be inaccurate")
#preferences on missings values
if (!is.null(na.act)) {
subset <- !(1:length(listw$neighbours) %in% na.act)
listw <- subset(listw, subset, zero.policy = zero.policy)
}
#check that W ius an object of class listw
if(!inherits(listw,"listw"))
stop("The weights matrix is not a listw object")
#generates x and y
y<-model.extract(mf,"response")
x<-model.matrix(mt,mf)
#checks on teh dimensions of x and y
if (length(y)!=nrow(x))
stop("x and y have different length")
#check on the dimensions of x and W
if (nrow(x) != length(listw$neighbours))
stop("Input data and weights have different dimension")
#check that X and y does not have missing values
if (any(is.na(y)))
stop("NAs in dependent variable")
if (any(is.na(x)))
stop("NAs in independent variable")
#fix the dimensions of the problem
n<-nrow(x)
k<-ncol(x)
xcolnames<-colnames(x)
K<-ifelse(xcolnames[1] == "(Intercept)" || all(x[ ,1]==1), 2, 1)
if(sarar){
wy<-lag.listw(listw,y, zero.policy=zero.policy)
wy<-array(wy,c(length(y),1))
colnames(wy)<-("Wy")
if (any(is.na(wy)))
stop("NAs in spatially lagged dependent variable")
if (k > 1) {
WX <- matrix(nrow = n, ncol = (k - (K - 1)))
WWX <- matrix(nrow = n, ncol = (k - (K - 1)))
for (i in K:k) {
wx <- lag.listw(listw, x[, i], zero.policy = zero.policy)
wwx<- lag.listw(listw, wx, zero.policy = zero.policy)
if (any(is.na(wx)))
stop("NAs in lagged independent variable")
WX[, (i - (K - 1))] <- wx
WWX[, (i - (K - 1))] <- wwx
}
}
instr<-cbind(WX,WWX)
##spatial two stage least square of the initial model
firststep<-tsls(y=y,yend=wy, X=x, Zinst = instr)
ubase<-residuals(firststep)
#GM step 1b
int1<-Ggfastfast(listw=listw,ubase,n, zero.policy = zero.policy)
##One could start from this initial value but if uses optimize, the initial values are not an issue
if (initial.value=="SAR"){
Wubase<-lag.listw(listw,ubase, zero.policy = zero.policy)
pars<-coefficients(lm(ubase~Wubase-1))
}
else pars<-initial.value
optres1<-nlminb(pars, arg, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1)
param<-optres1$par
# print(param)
reg<-x
u<-ubase
toinst<-wy
fi<-fifour(x,listw,instr,ubase,wy,param,n,inverse,eps, zero.policy= zero.policy)
# print(fi$nlw)
GMMfeas3<-nlminb(param, arg1, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1, VC=fi$nlw)
rhotilde <- GMMfeas3$par
# print(rhotilde)
fi1<-fifour(x,listw,instr,ubase,wy,rhotilde,n,inverse=inverse,eps=eps, zero.policy = zero.policy)
Om<-Omega(n,gamma=fi1$res1,H=fi1$instr,param=rhotilde,G=int1$bigG,FIinv=fi1$nlw,a=fi1$a,FI=fi1$nl,P=fi1$P,gammas=fi1$gammas,Ws=fi1$Ws)
yt <- y - rhotilde * wy
xt <- x - rhotilde * lag.listw(listw,x, zero.policy = zero.policy)
wyt<-wy - rhotilde * lag.listw(listw,wy, zero.policy = zero.policy)
colnames(xt)<-xcolnames
colnames(wyt)<-c('Wyt')
secstepb<-tsls(y=yt, yend=wyt, X=xt, Zinst = instr, reg=x, end=wy, yor=y,modified=TRUE)
utildeb<-secstepb$residuals
int1b<-Ggfastfast(listw, utildeb,n, zero.policy = zero.policy)
psippb<-fistslsfast(reg=xt, Ws=fi$Ws, instr=instr, resid=utildeb, toinst=wyt, param=rhotilde,solo=x,end=wy)
pars1<-rhotilde
GMMstsls1b<-nlminb(pars1, arg1, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1b, VC=psippb$nlw)
rhohat<-GMMstsls1b$par
#xt <- x - rhohat * lag.listw(listw,x)
psifinal<-fistslsfast(reg=xt, Ws=fi$Ws, instr=instr, resid=utildeb, toinst=wyt, param=rhohat,solo=x,end=wy)
Omfinal<- Omegabis(gammas=psifinal$res1, Hs=psifinal$instr, param=rhohat, G=int1b$bigG, FIinv=psifinal$nlw, as=psifinal$a, n, FI=psifinal$nl, P=psifinal$P)
coef<-(secstepb$coefficients[2:(ncol(xt)+1)])
coeff<-as.matrix(c(coef, secstepb$coefficients[1], rhohat ))
rownames(coeff)<-c(xcolnames, 'lambda', 'rho')
s2<-crossprod(utildeb)/(n-k)
model.data<-data.frame(cbind(y,x[,-1]))
method<-"gs2slshac"
k<-nrow(coeff)
R<-matrix(0,1,k)
R[,((k-1):k)]<-1
Rbeta<-R%*%coeff
Rvar<-R%*%Omfinal$VarCov%*%t(R)
stat<-as.numeric(t(Rbeta)%*% Rbeta/Rvar)
pval <- pchisq(stat,df=1,lower.tail=FALSE)
W<-list(stat=stat,pval=pval)
results<-list(coefficients=coeff,var=Omfinal$VarCov, s2=s2, call=cl, residuals=as.numeric(utildeb), model=model.data,method=method,W=W,yhat=secstepb$yhat )
}
else{
##this is an error model estimated by ols
firststep<-lm(y~x-1)
ubase<-residuals(firststep)
int1<-Ggfastfast(listw=listw,ubase,n, zero.policy = zero.policy)
##One could start from this initial value but if uses optimize, the initial values are not an issue
if (initial.value=="SAR"){
Wubase<-lag.listw(listw,ubase, zero.policy = zero.policy)
pars<-coefficients(lm(ubase~Wubase-1))
}
else pars<-initial.value
optres1<-nlminb(pars, arg, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1)
param<-optres1$par
reg<-x
u<-ubase
fi<-fierror(x,listw, ubase, param,n,inverse=inverse,eps=eps)
GMMfeas3<-nlminb(param, arg1, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1, VC=fi$nlw)
rhotilde<- GMMfeas3$par
# print(rhotilde)
fi1<-fierror(x,listw,ubase,rhotilde,n,inverse=inverse,eps=eps)
Om<-Omega(n,gamma=fi1$res1, H=x, param=rhotilde,G=int1$bigG,FIinv=fi1$nlw,a=fi1$a,FI=fi1$nl,P=fi1$P,gammas=fi1$gammas,Ws=fi1$Ws) ## in an error model H=X
#####################
yt <- y - rhotilde * lag.listw(listw,y, zero.policy=zero.policy)
xt <- x - rhotilde * lag.listw(listw,x, zero.policy=zero.policy)
colnames(xt)<-xcolnames
secstepb<-lm(yt~xt-1)
utildeb<-y - x %*% as.matrix(coefficients(secstepb))
int1b<-Ggfastfast(listw, utildeb,n, zero.policy = zero.policy)
psippb<-fistslserror(reg=x, Ws=fi$Ws, resid=utildeb, param=rhotilde, solo = cbind(x,lag.listw(listw,x[,-1], zero.policy=zero.policy)))
pars1<-rhotilde
GMMstsls1b<-nlminb(pars1, arg1, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1b, VC=psippb$nlw)
rhohat<-GMMstsls1b$par
psifinal<-fistslserror(reg=x, Ws=fi$Ws, resid=utildeb, param=rhohat,solo=cbind(x,lag.listw(listw,x[,-1], zero.policy=zero.policy)))
Omfinal<- Omegabis(gammas=psifinal$res1, Hs = cbind(x,lag.listw(listw,x[,-1], zero.policy=zero.policy)) , param=rhohat, G=int1b$bigG, FIinv=psifinal$nlw, as=psifinal$a, n, FI=psifinal$nl, P=psifinal$P)
coef<-coefficients(secstepb)
coeff<-as.matrix(c(coef, rhohat ))
rownames(coeff)<-c(xcolnames, 'rho')
s2<-crossprod(utildeb)/(n-k)
model.data<-data.frame(cbind(y,x[,-1]))
method<-"gs2slshac"
results<-list(coefficients=coeff,var=Omfinal$VarCov, s2=s2, call=cl, residuals=as.numeric(utildeb), model=model.data,method=method,W=NULL,yhat=fitted(secstepb) )
}
class(results)<-c("sphet", "gstsls")
return(results)
}
impacts.gstsls <- function(obj, ..., tr=NULL, R=NULL, listw=NULL,
evalues=NULL, tol=1e-6, empirical=FALSE, Q=NULL) {
if (!is.null(obj$listw_style) && obj$listw_style != "W")
stop("Only row-standardised weights supported")
coefs <- drop(obj$coefficients)
p2 <- length(coefs)
rho <- coefs[(p2-1)]
beta <- coefs[1:(p2-2)]
lambda <- coefs[p2]
# rho is lag coef., lambda is error coef (reversed from function use)
icept <- grep("(Intercept)", names(beta))
iicept <- length(icept) > 0
if (iicept) {
P <- matrix(beta[-icept], ncol=1)
bnames <- names(beta[-icept])
} else {
P <- matrix(beta, ncol=1)
bnames <- names(beta)
}
p <- length(beta)
n <- length(obj$residuals)
mu <- c(beta, rho, lambda)
Sigma <- obj$var
irho <- p2-1
drop2beta <- c(p2-1, p2)
if (!requireNamespace("spatialreg", quietly=TRUE))
stop("install spatialreg")
res <- spatialreg::intImpacts(rho=rho, beta=beta, P=P, n=n, mu=mu,
Sigma=Sigma, irho=irho, drop2beta=drop2beta, bnames=bnames,
interval=NULL, type="lag", tr=tr, R=R, listw=listw, evalues=evalues,
tol=tol, empirical=empirical, Q=Q, icept=icept, iicept=iicept, p=p)
attr(res, "iClass") <- class(obj)
res
}
Try the sphet package in your browser
Any scripts or data that you put into this service are public.
sphet documentation built on Jan. 6, 2022, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions impacts.gstsls gstslshet
Documented in gstslshet impacts.gstsls
#' GM estimation of a Cliff-Ord type model with Heteroskedastic Innovations
#' @aliases gstslshet
#' @name gstslshet
#'
#' @description Multi step GM/IV estimation of a linear Cliff and Ord -type of model
#' of the form:
#'
#' \deqn{y=\lambda W y + X \beta + u }
#' \deqn{u=\rho W u + e} with \deqn{e ~ N(0,\sigma^2_i) }
#'
#' The model allows for spatial lag in the dependent variable
#' and disturbances. The innovations in the disturbance process are assumed
#' heteroskedastic of an unknown form.
#'
#' @usage gstslshet(formula, data = list(), listw, na.action = na.fail,
#' zero.policy = NULL, initial.value = 0.2, abs.tol = 1e-20,
#' rel.tol = 1e-10, eps = 1e-5, inverse = T, sarar = T)
#'
#' @param formula a description of the model to be fit
#' @param data an object of class \link{data.frame}. An optional data frame containing the variables
#' in the model.
#' @param listw an object of class \code{listw} created for example by \code{nb2listw}
#' @param na.action a function which indicates what should happen when the data contains missing values.
#' See \link{lm} for details.
#' @param zero.policy See \code{lagsarlm} for details
#' @param initial.value The initial value for \eqn{\rho}. It can be either numeric (default is 0.2) or
#' set to \code{'SAR'}, in which case the optimization will start from the estimated coefficient of a regression of the 2SLS
#' residuals over their spatial lag (i.e. a spatial AR model)
#' @param abs.tol Absolute tolerance. See \link{nlminb} for details.
#' @param rel.tol Relative tolerance. See \link{nlminb} for details.
#' @param eps Tolerance level for the approximation. See Details.
#' @param inverse \code{TRUE}. If \code{FALSE}, an appoximated inverse is calculated. See Details.
#' @param sarar \code{TRUE}. If \code{FALSE}, a spatial error model is estimated.
#'
#' @details The procedure consists of two steps alternating GM and IV estimators. Each step consists of sub-steps.
#' In step one \eqn{\delta = [\beta',\lambda]'} is estimated by 2SLS. The 2SLS residuals are first employed
#' to obtain an initial (consistent but not efficient) GM estimator of \eqn{\rho} and then a consistent and efficient
#' estimator (involving the variance-covariance matrix of the limiting distribution of the normalized sample moments).
#' In step two, the spatial Cochrane-Orcutt transformed model is estimated by 2SLS. This corresponds to a GS2SLS procedure.
#' The GS2SLS residuals are used to obtain a consistent and efficient GM estimator for \eqn{\rho}.
#'
#' The initial value for the optimization in step 1b is taken to be \code{initial.value}. The initial value in step 1c is the
#' optimal parameter of step 1b. Finally, the initial value for the optimization of step 2b is the optimal parameter of step 1c.
#'
#' Internally, the object of class \code{listw} is transformed into a \link{Matrix}
#' using the function \link{listw2dgCMatrix}.
#'
#'
#' The expression of the estimated variance covariance matrix of the limiting
#' distribution of the normalized sample moments based on 2SLS residuals
#' involves the inversion of \eqn{I-\rho W'}.
#' When \code{inverse} is \code{FALSE}, the inverse is calculated using the approximation
#' \eqn{I +\rho W' + \rho^2 W'^2 + ...+ \rho^n W'^n}.
#' The powers considered depend on a condition.
#' The function will keep adding terms until the absolute value of the \code{sum} of all elements
#' of the matrix \eqn{\rho^i W^i} is greater than a fixed \eqn{\epsilon} (\code{eps}). By default \code{eps}
#' is set to 1e-5.
#'
#'
#' @return A list object of class \code{sphet}
#' \item{coefficients}{Generalized Spatial two stage least squares coefficient estimates of \eqn{\delta} and GM estimator for \eqn{\rho}. }
#' \item{var}{variance-covariance matrix of the estimated coefficients}
#' \item{s2}{GS2SLS residuals variance}
#' \item{residuals}{GS2SLS residuals}
#' \item{yhat}{difference between GS2SLS residuals and response variable}
#' \item{call}{the call used to create this object}
#' \item{model}{the model matrix of data}
#' \item{method}{\code{'gs2slshac'}}
#' \item{W}{Wald test for both \eqn{\rho} and \eqn{\lambda} are zero}
#'
#' @references
#' Arraiz, I. and Drukker, M.D. and Kelejian, H.H. and Prucha, I.R. (2007)
#' A spatial Cliff-Ord-type Model with Heteroskedastic Innovations: Small and Large Sample Results,
#' \emph{Department of Economics, University of Maryland}'
#'
#' Kelejian, H.H. and Prucha, I.R. (2007)
#'Specification and Estimation of Spatial Autoregressive Models with Autoregressive and Heteroskedastic Disturbances,
#' \emph{Journal of Econometrics}, forthcoming.
#'
#' Kelejian, H.H. and Prucha, I.R. (1999)
#' A Generalized Moments Estimator for the Autoregressive Parameter in a Spatial Model,
#' \emph{International Economic Review}, \bold{40}, pages 509--533.
#'
#' Kelejian, H.H. and Prucha, I.R. (1998)
#'A Generalized Spatial Two Stage Least Square Procedure for Estimating a Spatial Autoregressive
#' Model with Autoregressive Disturbances,
#' \emph{Journal of Real Estate Finance and Economics}, \bold{17}, pages 99--121.
#'
#' @seealso \code{\link{stslshac}}
#'
#' @author Gianfranco Piras \email{gpiras@mac.com}
#'
#' @examples
#' data(columbus, package = "spdep")
#' listw <- spdep::nb2listw(col.gal.nb)
#' res <- gstslshet(CRIME ~ HOVAL + INC, data = columbus, listw = listw)
#' summary(res)
#'
#' @keywords models
#' @export
gstslshet<-function(formula, data=list(),listw, na.action=na.fail, zero.policy=NULL, initial.value=0.2, abs.tol=1e-20, rel.tol=1e-10, eps=1e-5, inverse=T, sarar=T){
##functions that need to be sourced
#source("twostagels.R")
#source("utilities.R")
#source("listw2dgCMatrix.R")
#source("Omega.R")
if (is.null(zero.policy))
zero.policy <- get.ZeroPolicyOption()
stopifnot(is.logical(zero.policy))
#extract model objects
mt<-terms(formula,data=data)
mf<-lm(formula, data, na.action=na.action, method="model.frame")
na.act<-attr(mf,'na.action')
# record call
cl<- match.call()
#if(!inverse) warning("Approximated inverse could be inaccurate")
#preferences on missings values
if (!is.null(na.act)) {
subset <- !(1:length(listw$neighbours) %in% na.act)
listw <- subset(listw, subset, zero.policy = zero.policy)
}
#check that W ius an object of class listw
if(!inherits(listw,"listw"))
stop("The weights matrix is not a listw object")
#generates x and y
y<-model.extract(mf,"response")
x<-model.matrix(mt,mf)
#checks on teh dimensions of x and y
if (length(y)!=nrow(x))
stop("x and y have different length")
#check on the dimensions of x and W
if (nrow(x) != length(listw$neighbours))
stop("Input data and weights have different dimension")
#check that X and y does not have missing values
if (any(is.na(y)))
stop("NAs in dependent variable")
if (any(is.na(x)))
stop("NAs in independent variable")
#fix the dimensions of the problem
n<-nrow(x)
k<-ncol(x)
xcolnames<-colnames(x)
K<-ifelse(xcolnames[1] == "(Intercept)" || all(x[ ,1]==1), 2, 1)
if(sarar){
wy<-lag.listw(listw,y, zero.policy=zero.policy)
wy<-array(wy,c(length(y),1))
colnames(wy)<-("Wy")
if (any(is.na(wy)))
stop("NAs in spatially lagged dependent variable")
if (k > 1) {
WX <- matrix(nrow = n, ncol = (k - (K - 1)))
WWX <- matrix(nrow = n, ncol = (k - (K - 1)))
for (i in K:k) {
wx <- lag.listw(listw, x[, i], zero.policy = zero.policy)
wwx<- lag.listw(listw, wx, zero.policy = zero.policy)
if (any(is.na(wx)))
stop("NAs in lagged independent variable")
WX[, (i - (K - 1))] <- wx
WWX[, (i - (K - 1))] <- wwx
}
}
instr<-cbind(WX,WWX)
##spatial two stage least square of the initial model
firststep<-tsls(y=y,yend=wy, X=x, Zinst = instr)
ubase<-residuals(firststep)
#GM step 1b
int1<-Ggfastfast(listw=listw,ubase,n, zero.policy = zero.policy)
##One could start from this initial value but if uses optimize, the initial values are not an issue
if (initial.value=="SAR"){
Wubase<-lag.listw(listw,ubase, zero.policy = zero.policy)
pars<-coefficients(lm(ubase~Wubase-1))
}
else pars<-initial.value
optres1<-nlminb(pars, arg, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1)
param<-optres1$par
# print(param)
reg<-x
u<-ubase
toinst<-wy
fi<-fifour(x,listw,instr,ubase,wy,param,n,inverse,eps, zero.policy= zero.policy)
# print(fi$nlw)
GMMfeas3<-nlminb(param, arg1, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1, VC=fi$nlw)
rhotilde <- GMMfeas3$par
# print(rhotilde)
fi1<-fifour(x,listw,instr,ubase,wy,rhotilde,n,inverse=inverse,eps=eps, zero.policy = zero.policy)
Om<-Omega(n,gamma=fi1$res1,H=fi1$instr,param=rhotilde,G=int1$bigG,FIinv=fi1$nlw,a=fi1$a,FI=fi1$nl,P=fi1$P,gammas=fi1$gammas,Ws=fi1$Ws)
yt <- y - rhotilde * wy
xt <- x - rhotilde * lag.listw(listw,x, zero.policy = zero.policy)
wyt<-wy - rhotilde * lag.listw(listw,wy, zero.policy = zero.policy)
colnames(xt)<-xcolnames
colnames(wyt)<-c('Wyt')
secstepb<-tsls(y=yt, yend=wyt, X=xt, Zinst = instr, reg=x, end=wy, yor=y,modified=TRUE)
utildeb<-secstepb$residuals
int1b<-Ggfastfast(listw, utildeb,n, zero.policy = zero.policy)
psippb<-fistslsfast(reg=xt, Ws=fi$Ws, instr=instr, resid=utildeb, toinst=wyt, param=rhotilde,solo=x,end=wy)
pars1<-rhotilde
GMMstsls1b<-nlminb(pars1, arg1, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1b, VC=psippb$nlw)
rhohat<-GMMstsls1b$par
#xt <- x - rhohat * lag.listw(listw,x)
psifinal<-fistslsfast(reg=xt, Ws=fi$Ws, instr=instr, resid=utildeb, toinst=wyt, param=rhohat,solo=x,end=wy)
Omfinal<- Omegabis(gammas=psifinal$res1, Hs=psifinal$instr, param=rhohat, G=int1b$bigG, FIinv=psifinal$nlw, as=psifinal$a, n, FI=psifinal$nl, P=psifinal$P)
coef<-(secstepb$coefficients[2:(ncol(xt)+1)])
coeff<-as.matrix(c(coef, secstepb$coefficients[1], rhohat ))
rownames(coeff)<-c(xcolnames, 'lambda', 'rho')
s2<-crossprod(utildeb)/(n-k)
model.data<-data.frame(cbind(y,x[,-1]))
method<-"gs2slshac"
k<-nrow(coeff)
R<-matrix(0,1,k)
R[,((k-1):k)]<-1
Rbeta<-R%*%coeff
Rvar<-R%*%Omfinal$VarCov%*%t(R)
stat<-as.numeric(t(Rbeta)%*% Rbeta/Rvar)
pval <- pchisq(stat,df=1,lower.tail=FALSE)
W<-list(stat=stat,pval=pval)
results<-list(coefficients=coeff,var=Omfinal$VarCov, s2=s2, call=cl, residuals=as.numeric(utildeb), model=model.data,method=method,W=W,yhat=secstepb$yhat )
}
else{
##this is an error model estimated by ols
firststep<-lm(y~x-1)
ubase<-residuals(firststep)
int1<-Ggfastfast(listw=listw,ubase,n, zero.policy = zero.policy)
##One could start from this initial value but if uses optimize, the initial values are not an issue
if (initial.value=="SAR"){
Wubase<-lag.listw(listw,ubase, zero.policy = zero.policy)
pars<-coefficients(lm(ubase~Wubase-1))
}
else pars<-initial.value
optres1<-nlminb(pars, arg, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1)
param<-optres1$par
reg<-x
u<-ubase
fi<-fierror(x,listw, ubase, param,n,inverse=inverse,eps=eps)
GMMfeas3<-nlminb(param, arg1, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1, VC=fi$nlw)
rhotilde<- GMMfeas3$par
# print(rhotilde)
fi1<-fierror(x,listw,ubase,rhotilde,n,inverse=inverse,eps=eps)
Om<-Omega(n,gamma=fi1$res1, H=x, param=rhotilde,G=int1$bigG,FIinv=fi1$nlw,a=fi1$a,FI=fi1$nl,P=fi1$P,gammas=fi1$gammas,Ws=fi1$Ws) ## in an error model H=X
#####################
yt <- y - rhotilde * lag.listw(listw,y, zero.policy=zero.policy)
xt <- x - rhotilde * lag.listw(listw,x, zero.policy=zero.policy)
colnames(xt)<-xcolnames
secstepb<-lm(yt~xt-1)
utildeb<-y - x %*% as.matrix(coefficients(secstepb))
int1b<-Ggfastfast(listw, utildeb,n, zero.policy = zero.policy)
psippb<-fistslserror(reg=x, Ws=fi$Ws, resid=utildeb, param=rhotilde, solo = cbind(x,lag.listw(listw,x[,-1], zero.policy=zero.policy)))
pars1<-rhotilde
GMMstsls1b<-nlminb(pars1, arg1, lower= -0.9 + .Machine$double.eps , upper= 0.9 - .Machine$double.eps,control=list(abs.tol=abs.tol,rel.tol=rel.tol),v=int1b, VC=psippb$nlw)
rhohat<-GMMstsls1b$par
psifinal<-fistslserror(reg=x, Ws=fi$Ws, resid=utildeb, param=rhohat,solo=cbind(x,lag.listw(listw,x[,-1], zero.policy=zero.policy)))
Omfinal<- Omegabis(gammas=psifinal$res1, Hs = cbind(x,lag.listw(listw,x[,-1], zero.policy=zero.policy)) , param=rhohat, G=int1b$bigG, FIinv=psifinal$nlw, as=psifinal$a, n, FI=psifinal$nl, P=psifinal$P)
coef<-coefficients(secstepb)
coeff<-as.matrix(c(coef, rhohat ))
rownames(coeff)<-c(xcolnames, 'rho')
s2<-crossprod(utildeb)/(n-k)
model.data<-data.frame(cbind(y,x[,-1]))
method<-"gs2slshac"
results<-list(coefficients=coeff,var=Omfinal$VarCov, s2=s2, call=cl, residuals=as.numeric(utildeb), model=model.data,method=method,W=NULL,yhat=fitted(secstepb) )
}
class(results)<-c("sphet", "gstsls")
return(results)
}
impacts.gstsls <- function(obj, ..., tr=NULL, R=NULL, listw=NULL,
evalues=NULL, tol=1e-6, empirical=FALSE, Q=NULL) {
if (!is.null(obj$listw_style) && obj$listw_style != "W")
stop("Only row-standardised weights supported")
coefs <- drop(obj$coefficients)
p2 <- length(coefs)
rho <- coefs[(p2-1)]
beta <- coefs[1:(p2-2)]
lambda <- coefs[p2]
# rho is lag coef., lambda is error coef (reversed from function use)
icept <- grep("(Intercept)", names(beta))
iicept <- length(icept) > 0
if (iicept) {
P <- matrix(beta[-icept], ncol=1)
bnames <- names(beta[-icept])
} else {
P <- matrix(beta, ncol=1)
bnames <- names(beta)
}
p <- length(beta)
n <- length(obj$residuals)
mu <- c(beta, rho, lambda)
Sigma <- obj$var
irho <- p2-1
drop2beta <- c(p2-1, p2)
if (!requireNamespace("spatialreg", quietly=TRUE))
stop("install spatialreg")
res <- spatialreg::intImpacts(rho=rho, beta=beta, P=P, n=n, mu=mu,
Sigma=Sigma, irho=irho, drop2beta=drop2beta, bnames=bnames,
interval=NULL, type="lag", tr=tr, R=R, listw=listw, evalues=evalues,
tol=tol, empirical=empirical, Q=Q, icept=icept, iicept=iicept, p=p)
attr(res, "iClass") <- class(obj)
res
}
Try the sphet package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.