Nothing
KSS.default <- function(formula,
additive.effects = c("none", "individual", "time", "twoways"),
consult.dim.crit = FALSE,
d.max = NULL,
sig2.hat = NULL,
factor.dim = NULL,
level = 0.01,
spar = NULL,
CV = FALSE,
convergence = 1e-6,
restrict.mode = c("restrict.factors","restrict.loadings"),
...)
{
##===================================================================================
if(!class(formula)=="formula"){
stop("\n Argument >>formula<< needs a formula-object like y~x1+... where the elements are matrices.")
}
if(!any(additive.effects==c("none", "individual", "time", "twoways"))){
stop("\n Argument >>effect<< must be one of: \n none, individual, time, twoways")
}
if(!is.numeric(level)){
stop("\n Argument >>alpha<< has to be numeric.")
}
if(!is.null(factor.dim) & !is.numeric(factor.dim)){
stop("\n Argument >>factor.dim<< has to be numeric.")
}
if(!is.null(spar) & !is.numeric(spar)){
stop("\n Argument >>spar<< has to be numeric.")
}
if(!is.logical(CV)){
stop("\n Argument >>CV<< has to be TRUE or FALSE.")
}
if(!is.numeric(convergence)){
stop("\n Argument >>convergence<< has to be numeric.")
}
if(!is.null(sig2.hat) &!is.numeric(sig2.hat)){
stop("\n Argument >>sig2.hat<< has to be numeric.")
}
if(!any(restrict.mode==c("restrict.factors","restrict.loadings"))){
stop("\n Argument >>restrict.mode<< must be either: restrict.factors or: restrict.loadings")
}
##====================================================================================
## check "effect" and "dim.crit"
effect <- match.arg(additive.effects)
## extract data from formula
names <- names(model.frame(formula))
PF.obj <- FUN.Pformula(formula = formula, effect = effect)
N <- ncol(PF.obj[[1]]$ODM)
T <- nrow(PF.obj[[1]]$ODM)
P <- length(PF.obj)-1
dat.dim <- c(T, N, P)
is.intercept <- PF.obj[[1]]$I
## *OR*iginal *dat*a
ORdat <- sapply(1:(P+1), function(i) PF.obj[[i]]$ODM)
## *TR*ansformed *dat*a
TRdat <- sapply(1:(P+1), function(i) PF.obj[[i]]$TDM)
dat.matrix <- sapply(1:(P+1), function(i) PF.obj[[i]]$TDV)
Or.Y <- ORdat[, 1, drop = FALSE] # (TN x 1)
Or.X <- ORdat[, 2:(P+1), drop = FALSE] # (TN x P)
TR.Y <- TRdat[, 1, drop = FALSE] # (TN x 1)
TR.X <- TRdat[, 2:(P+1), drop = FALSE] # (TN x P)
TR.Y.mat <- matrix(TR.Y, T, N) # (T x N)
TR.X.mat <- matrix(TR.X, T, (N*P)) # (T x NP)
#####################################################################################################
## GCV Procedure to determine the smoothing parameter
iterateGCV <- FUN.iterate.GCV(TR.Y.mat, TR.X.mat, N, T, P)
spar.GCV <- iterateGCV$PCA$spar
if(is.null(spar)){
spar.low <- spar.GCV*0.75
}else{
spar.low <- spar
}
if(CV){
spar.interval.max <- spar.GCV
spar.CV <- KSS.CV(kappa.interv=c(.Machine$double.eps, spar.interval.max),
Y=TR.Y, X=TR.X, N=N, T=T, P=P, spar.dim.fit=spar.low, tol=convergence)$minimum
cat("\n CV-Optimization converged.\n")
spar.low <- spar.CV
}
#####################################################################################################
t.seq <- seq(0, 1, length.out=T)
## Determination of the 'norder' argument of function 'smooth.Pspline()'
if(T>=5){
norder <- 2
}
if(T==4){
norder <- 1
}
if(T<4){
stop("Time-Dimension T of the panel data has to be at least T=4.")
}
TR.Y.mat.smth <- smooth.Pspline(x = t.seq, y = TR.Y.mat, norder = norder, spar = spar.low)$ysmth #(T x N)
TR.X.mat.smth <- smooth.Pspline(x = t.seq, y = TR.X.mat, norder = norder, spar = spar.low)$ysmth #(T x NP)
TR.X.mat.smth2 <- smooth.Pspline(x = t.seq, y = TR.X.mat.smth, norder = norder, spar = spar.low)$ysmth #(T x NP)
## calculate beta coefficents
TR.Y.smth <- matrix(TR.Y.mat.smth, nrow= (N*T), ncol = 1) # (TN x 1)
TR.X.smth <- matrix(TR.X.mat.smth, nrow= (N*T), ncol = P) # (TN x P)
TR.X.smth2 <- matrix(TR.X.mat.smth2, nrow= (N*T), ncol = P) # (TN x P)
t.TR.X.TR.X <- crossprod(TR.X) # (PxP)
t.TR.X.TR.X.smth <- crossprod(TR.X, TR.X.smth) # (PxP)
t.TR.X.TR.X.smth2<- crossprod(TR.X, TR.X.smth2) # (PxP)
t.TR.X.TR.Y <- crossprod(TR.X, TR.Y) # (Px1)
t.TR.X.TR.Y.smth <- crossprod(TR.X, TR.Y.smth) # (Px1)
bloc1 <- t.TR.X.TR.X - t.TR.X.TR.X.smth # (PxP)
bloc2 <- t.TR.X.TR.Y - t.TR.X.TR.Y.smth # (Px1)
## common-Slope.Coefficients:
inv.bloc1 <- solve(bloc1)
com.slops.0 <- inv.bloc1%*%bloc2 # (Px1)
## calculate first step residuals and estimate dimension of factor-structure
Residu.mat <- matrix((TR.Y - (TR.X %*% com.slops.0)), T, N)
## functional pca
fpca.fit.obj <- fpca.fit(Residu.mat, spar=spar.low)
## Estimation of Dimension
dim.criterion <- c("PC1", "PC2", "PC3", "BIC3","IC1", "IC2", "IC3",
"IPC1","IPC2", "IPC3",
"ABC.IC1", "ABC.IC2",
"KSS.C",
"ED",
"ER", "GR")
Opt.dim.Output <- as.matrix(sapply(dim.criterion, function(dim.criterion){
EstDim(dim.criterion, Obj=Residu.mat, d.max=d.max, factor.dim=factor.dim,
sig2.hat=sig2.hat, level=level)[2]}))
## User-Interface---------------------------------------------------------------------------
Opt.dim.Output.BaiNg <- c(as.numeric(Opt.dim.Output[1:7,1]))
names(Opt.dim.Output.BaiNg) <- c("PC1","PC2","PC3","BIC3","IC1","IC2","IC3")
Opt.dim.Output.Bai <- c(as.numeric(Opt.dim.Output[8:10,1]))
names(Opt.dim.Output.Bai) <- c("IPC1","IPC2","IPC3")
Opt.dim.Output.Alessi <- c(as.numeric(Opt.dim.Output[11:12,1]))
names(Opt.dim.Output.Alessi) <- c("ABC.IC1", "ABC.IC2")
Opt.dim.Output.KSS <- c(as.numeric(Opt.dim.Output[13,1]))
names(Opt.dim.Output.KSS) <- c(" KSS.C")
Opt.dim.Output.Onatski <- c(as.numeric(Opt.dim.Output[11,1]))
names(Opt.dim.Output.Onatski) <- c(" ED")
Opt.dim.Output.RH <- c(as.numeric(Opt.dim.Output[12:13,1]))
names(Opt.dim.Output.RH) <- c(" ER","GR")
if(is.null(factor.dim) && consult.dim.crit){
cat("-----------------------------------------------------------\n")
cat("Results of Dimension-Estimations");
cat("\n\n-Bai and Ng (2002):\n"); print(Opt.dim.Output.BaiNg, quote = FALSE, na.print="");
cat("\n\n-Bai (2004):\n"); print(Opt.dim.Output.Bai, quote = FALSE, na.print="");
cat("\n\n-Alessi et al. (2010):\n"); print(Opt.dim.Output.Alessi, quote = FALSE, na.print="");
cat("\n-Kneip et al. (2012):\n"); print(Opt.dim.Output.KSS, quote = FALSE, na.print="");
cat("\n-Onatski (2009):\n"); print(Opt.dim.Output.Onatski, quote = FALSE, na.print="");
cat("\n-Ahn and Horenstein (2013):\n"); print(Opt.dim.Output.RH, quote = FALSE, na.print="");
otp <- as.numeric(c(Opt.dim.Output.BaiNg,Opt.dim.Output.KSS,Opt.dim.Output.Onatski,Opt.dim.Output.RH))
cat("\n")
cat("-----------------------------------------------------------\n")
## Auxiliary Fct's
myASK <- function(){
cat("Please, choose one of the proposed integers: ")
readLines(con = stdin(), n = 1)
}
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5){
if(is.na(x)){
result <- FALSE
}else{
if(is.numeric(x) & abs(x - round(x)) < tol){
result <- min(as.numeric(Opt.dim.Output[,1]))<=x & x<= max(as.numeric(Opt.dim.Output[,1]))
}else{
result <- FALSE
}
}
return(result)
}
used.dim <- suppressWarnings(as.numeric(myASK()))
## Checking the Comand-Line input:
while(!is.wholenumber(used.dim)){
cat("Wrong input! You have to select a particular dimension by giving a numeric integer value between ", min(as.numeric(Opt.dim.Output[,1]))," and ", max(as.numeric(Opt.dim.Output[,1])),"! \n")
used.dim <- suppressWarnings(as.numeric(myASK()))
}
cat("Used dimension of unobs. factor structure is:", used.dim,"\n")
cat("-----------------------------------------------------------\n")
}
##------------------------------------------------------------------------------------------
if(!is.null(factor.dim)){
used.dim <- factor.dim
}
if(is.null(factor.dim) && !consult.dim.crit){
used.dim <- c(as.numeric(Opt.dim.Output[13,1]))# KSS.C
}
## now: 'used.dim' is specified
if(used.dim > 0){
factors <- fpca.fit.obj$factors[, 1:used.dim, drop= FALSE]
loadings <- fpca.fit.obj$loadings[, 1:used.dim, drop= FALSE]
factor.stract <- tcrossprod(factors, loadings)
## Eventually for later extensions: logical argument two.step in order
## to allow for a conditional estimation of the slope-parameters given
## the estimated dimension d.
two.step <- FALSE
if(two.step){
## re-estimate beta=========================================================
NEW.TR.Y.mat <- TR.Y.mat - factor.stract
NEW.TR.Y <- as.vector(NEW.TR.Y.mat)
beta <- qr.solve(TR.X, NEW.TR.Y)
beta <- matrix(beta, P,1)
##==========================================================================
}else{
beta <- com.slops.0
}
}
if(used.dim == 0){# no fact.-struct
factors <- NULL
loadings <- NULL
factor.stract <- matrix(0,T,N)
beta <- qr.solve(TR.X, TR.Y) # OLS
}
## Built up the return object *est*
FUN.add.eff.obj <- FUN.add.eff(PF.obj = PF.obj,
fpca.fit.obj = fpca.fit.obj,
beta.hat = beta)
## degrees.of.freedom =======================================================================
degrees.of.freedom <- (T*N - (N+T)*used.dim - P -
is.intercept -
N*(effect == "individual" | effect == "twoways") -
T*(effect == "time" | effect == "twoways"))
## re-estimation of sig2.hat (Paper KSS Section 3.4)==========================================
YOVc <- PF.obj[[1]]$TRm$OVc
XOVc <- sapply(2:(P+1), function(i)PF.obj[[i]]$TRm$OVc)
Or.Y_minus_YOVc <- Or.Y - YOVc
Or.X_minus_XOVc <- Or.X - matrix(rep(XOVc, each=(N*T)), N*T, P)
sig2.hat <- sum((TR.Y - TR.X %*% beta - matrix(factor.stract, N*T, 1))^2)/degrees.of.freedom
## estimation of Variances ==================================================================
if(used.dim > 0){
## estimation of beta-variance beta.V
if(!two.step){
beta.V <- sig2.hat * solve(bloc1) %*% (t.TR.X.TR.X + t.TR.X.TR.X.smth2 - 2*t.TR.X.TR.X.smth) %*% solve(bloc1)
}
if(two.step){
beta.V <- sig2.hat * solve(crossprod(TR.X))
}
## estimation of Intercept-variance
if(is.intercept){
Intercept.V <- (sig2.hat + matrix(colMeans(Or.X),1,P) %*% beta.V %*% t(matrix(colMeans(Or.X),1,P)))
}else{Intercept.V <- NULL}
}
if(used.dim == 0){
if(!is.intercept){
beta.V <- sig2.hat * solve(t.TR.X.TR.X)
Intercept.V <- NULL
}else{
beta.V <- sig2.hat * solve(crossprod(TR.X))
Intercept.V <- (sig2.hat + matrix(colMeans(Or.X),1,P) %*% beta.V %*% t(matrix(colMeans(Or.X),1,P)))
}
}
## Fitted Values =============================================================================
fitted.values <- matrix(c(rep(FUN.add.eff.obj$mu, T*N) + rep(FUN.add.eff.obj$beta.0, N) +
rep(FUN.add.eff.obj$tau, each=T) + Or.X %*% beta)
, T, N) + factor.stract
## Return ====================================================================================
est <- vector("list")
est$dat.matrix <- dat.matrix
est$formula <- formula
est$dat.dim <- dat.dim
est$slope.para <- beta
est$beta.V <- beta.V
est$names <- names #Names of: dependent variable and regressors
est$is.intercept <- is.intercept #Intercept: TRUE or FALSE
est$additive.effects <- effect #Additive-Effect-Type
est$Intercept <- FUN.add.eff.obj$mu # Intercept
est$Intercept.V <- Intercept.V
est$Add.Ind.Eff <- FUN.add.eff.obj$tau # Add. indiv. Effects
est$Add.Tim.Eff <- FUN.add.eff.obj$beta.0 # Add. time Effects (beta.0-function)
est$unob.factors <- factors
est$ind.loadings <- loadings
est$unob.fact.stru <- factor.stract
est$used.dim <- used.dim
est$optimal.dim <- Opt.dim.Output
est$fitted.values <- fitted.values
est$orig.Y <- matrix(Or.Y, T, N)
est$residuals <- est$orig.Y - est$fitted.values
est$sig2.hat <- sig2.hat
est$degrees.of.freedom <- degrees.of.freedom
est$call <- match.call()
class(est) <- "KSS"
##=============================================================================================
return(est)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.