# R/Eup.R In phtt: Panel Data Analysis with Heterogeneous Time Trends

#### Documented in Eup.default

```####################### Entire Updated Estimation #########################
# Input:
#	  1. dat.matrix	=  is a matrix where the first colomn containing
#		 the(NTx1) vector of Ys and the remaining colomns containing
#		the (NTxP) vector of Xs
# 	  2. dat.dim         dat.dim[1] is nr (nbr of rows) and dat.dim[2]
#		 isnc (nbr of colomns) of the panel matrix Y.
#	  3. dim.criterion	= c("PC1", "PC2", "PC3", "IC1", "IC2" , "IC3"
#		, "IPC1", "IPC2", "IPC3", "KSS.C1", "KSS.C2", "ED", "ER", "GR")
#	  4. factor.dim	= the number of factors if it is known (standard is
#		NULL)
#	  5. d.max 		= the maximum number of factos (an argument needed
#		for the dimension selction)
#	  6. sig2.hat = is an argument needed for the dimension selction
#       7. level	= is an argument needed for the dimension selction
#	  8. spar	= is an argument needed for the dimension selction
#	  9. double.iteration =  logical argument, if TRUE the function will
#		estimate the optimal dimension in an outer iteration
#		the convergence of all the paramters is obtained in a double
#		iteration. If FALSE the dimension will be estimated parallelly
#		with beta, lambda and F in order to reduce the number of
#		computation (desadvantege: convergence to a local a local
#		optimum).
#		This argument will be neglected if factor.dim is specified.
# Output:
#	  1. \$PCA	= svd.pca object calculated at the optimal dimension
#			(or given factor.dim)
#	  2. \$beta	= optimal dimension according to the dimensionalty
#			criterion
#	  3. \$opt.d = the slope estimator of the observed regressors
#			(beta.eup)
#	  4. \$nbr.iterations = number of iteration
##########################################################################

FUN.Eup <- function(dat.matrix, dat.dim
, double.iteration = double.iteration
, dim.criterion, factor.dim, d.max, sig2.hat
, start.beta, max.iteration, convergence){
#### data
y 	<- dat.matrix[, 1, drop = FALSE]
x 	<- dat.matrix[,-1, drop = FALSE]

nr 	<- dat.dim[1]
nc	<- dat.dim[2]
P	<- dat.dim[3]# or ncol(x)

#### if  d.max not given then d.max will be setted to sqrt(min(N, T))
if(is.null(d.max)) d.max <- round(sqrt(min(nr, nc)))

#### start value #####
beta.0 <- start.beta
if(is.null(beta.0))
{
d.max.star <- round(sqrt(min(nr, nc*P)))
ymats <- matrix(y, nr, nc)
xmats <- matrix(x, nr, (nc*P))
zmats <- cbind(ymats, xmats)
svdComp <- svd.pca(zmats, given.d = d.max.star)     # new
svdComp2 <- svd.pca(xmats, given.d = d.max.star)
Gstart <- svdComp\$L[, 0:d.max.star, drop = FALSE]
Gx     <- svdComp2\$L[, 0:d.max.star, drop = FALSE]

#Gstartnrnc <- matrix(t(Gstart), (nr*nc), d.max, byrow = TRUE)
Korr <- cor(Gstart, Gx)^2
delta <- diag((1- apply(Korr, 1, max)))
MG <- diag(1, nr) - Gstart%*%delta%*%t(Gstart)

XMG <- MG%*%xmats
YMG <- MG%*%ymats
trx <- matrix(XMG, (nr*nc), P)
try <- matrix(YMG, (nr*nc), 1)
beta.0 <- coef(lm(try ~ -1 + trx ))
}
#### calculate the inverse once in order to reduce computations of the
#### iterated slope estimator
FUN.ols.beta <- function(updated.y, x, inv.xx.x){
beta <- tcrossprod(inv.xx.x, t(updated.y))
}

xx 	   <- crossprod(x)
inv.xx   <- solve(xx)
inv.xx.x  <- inv.xx%*%t(x)

#### define given.d = factor.dim and set factor.dim= d.max if it is null
given.d <- factor.dim
if(is.null(factor.dim)) factor.dim <- d.max

#### Inner Iteration

inner.iteration <- function(y, x, inv.xx.x =inv.xx.x
, beta.0 = beta.0
, factor.dim = factor.dim, d.max=d.max
, sig2.hat= sig2.hat
, double.iteration = double.iteration
, max.iteration = max.iteration, convergence = convergence
, past.iterations = past.iterations, i=1){
# Iteration (0): initial cumputations
# w.0 = y-x%*%beta.0
w.0 <- y - tcrossprod(x, t(beta.0))

# W.0: write w.0 in a matrix form
W.0 <- matrix(w.0, nr, nc)

# PCA.0: PCA.0 computation, OptDim.0 and y.fitted
PCA.0    <- svd.pca(W.0, given.d=factor.dim)
if(!double.iteration && is.null(given.d)){
OptDim.0   <- EstDim(PCA.0
, dim.criterion = dim.criterion
, d.max = (factor.dim ), sig2.hat= sig2.hat) ### this differs from the version of phtt

opt.dim.0  <- OptDim.0[,2]
factor.dim <- opt.dim.0
y.fitted.0 <- tcrossprod(PCA.0\$L[, 0:opt.dim.0
, drop = FALSE])%*%W.0
}

else y.fitted.0 = PCA.0\$Q.fit

# Iteration (+1)
# y.updated.0: update y.updated.0 = y - fs.0
y.updated.1 <- y -  c(y.fitted.0)

# beta.1: OLS.1 computation for the computed fs.0 in interation 0
beta.1 <- FUN.ols.beta(y.updated.1, x, inv.xx.x)

# convergence condition
if(all( abs((beta.0 - beta.1)) < convergence)| i == max.iteration){
Result <- list(PCA=PCA.0, beta=beta.0
,factor.dim = factor.dim, Nbr.Iterations = i)
Result
}

else inner.iteration(y=y, x=x, inv.xx.x =inv.xx.x
,beta.0 = beta.1,factor.dim = factor.dim
,d.max=d.max,sig2.hat= sig2.hat
,double.iteration = double.iteration
,max.iteration = max.iteration, convergence = convergence
,past.iterations = past.iterations, (i+1))
}

####  outer and inner iterateion

entire.iteration <- function(y=y, x=x, inv.xx.x =inv.xx.x
, beta.0 = beta.0, factor.dim = factor.dim
, d.max=d.max, sig2.hat= sig2.hat
, double.iteration = double.iteration
, max.iteration = max.iteration
, convergence = convergence
, past.iterations = 0, l =1){

# suppose the algorithm will converge, if not 'converges' will turn up to FALSE
converges = TRUE
# first inner iteration
In.Iter.0 <- inner.iteration(y=y, x=x, inv.xx.x =inv.xx.x
, beta.0 = beta.0, factor.dim = factor.dim
, d.max=d.max, sig2.hat= sig2.hat
, double.iteration = double.iteration
, max.iteration = max.iteration, convergence = convergence
, past.iterations = past.iterations
, i =1)
pca.d0 	  <- In.Iter.0\$PCA
beta.d0 	  <- In.Iter.0\$beta
opt.d0	  <- In.Iter.0\$factor.dim
nbr.iteration <- In.Iter.0\$Nbr.Iterations + past.iterations

# if double.iteration is TRUE select new opt.d iteratively
if(double.iteration && is.null(given.d)){
# 1 new optimal dimension
opt.dim1 <- EstDim(pca.d0, dim.criterion = dim.criterion
, d.max = (opt.d0), sig2.hat= sig2.hat) ### this is the last version Eup as in the submitted paper and differs from the first version of phtt
opt.d1  <- opt.dim1[,2]
# convergence condition
if(opt.d1>=opt.d0){
if(In.Iter.0\$Nbr.Iterations == max.iteration){
converges = FALSE
#warning(paste("The maximal number of inner iterations is achieved ", sep=""), call. = FALSE)
}
Result  <- list(y=y, x=x, dat.dim= dat.dim
, PCA = pca.d0, beta= beta.d0
, opt.d =opt.d0
, nbr.iterations= nbr.iteration, converges = converges)
Result
}
else entire.iteration(y=y, x=x, inv.xx.x =inv.xx.x
, beta.0 = beta.d0, factor.dim = opt.d1
, d.max=d.max, sig2.hat= sig2.hat
, double.iteration = double.iteration
, past.iterations = nbr.iteration
, max.iteration = max.iteration, convergence = convergence
, l = (l+1))
}
else {
if(nbr.iteration == max.iteration){
converges = FALSE
#warning(paste("The maximal number of inner iterations is achieved ", sep=""), call. = FALSE)
}
Result <- list(y=y, x=x, dat.dim= dat.dim
, PCA = pca.d0, beta= beta.d0, opt.d =opt.d0
, nbr.iterations=nbr.iteration, converges = converges)
Result
}
}

###### entire iteration result

Result	 <- entire.iteration(y=y, x=x, inv.xx.x =inv.xx.x
, beta.0 = beta.0 , factor.dim = factor.dim
, d.max=d.max, sig2.hat= sig2.hat
, double.iteration = double.iteration
, max.iteration = max.iteration, convergence = convergence
, past.iterations = 0, l =1)

Result
}

############################### Eup.default ###############################
# Input:
#	  1. Formel
#	  2. additive.effects = c("none", "individual", "time", "twoways")
#	  3. dim.criterion	= c("PC1", "PC2", "PC3", "IC1", "IC2"
#			, "IC3", "IPC1", "IPC2", "IPC3", "KSS.C1"
#			, "KSS.C2", "ED", "ER", "GR")
#
#	  4. factor.dim	= the number of factors if it is known
#			(standard is NULL)
#	  5. d.max 		= the maximum number of factos
#			(an argument needed for the dimension selction)
#	  6. sig2.hat 	= is an argument needed for the dimension selction
#       7. level		= is an argument needed for the dimension selction
#	  8. spar		= is an argument needed for the dimension selction
#	  9. double.iteration =  logical argument, if TRUE the function will
#			estimate the optimal dimension in an outer iteration the
#			convergence of all the paramters is obtained in a double
#			iteration. If FALSE the dimension will be estimated
#			parallelly with beta, lambda and F in order to reduce
#			the number of computation (desadvantege: convergence to
#			a local a local optimum). This argument will be neglected
#			if factor.dim is specified.
# Output:
#	  1. \$PCA   		= svd.pca object calculated at the optimal
#					dimension (or given factor.dim)
#	  2. \$beta  		= the slope estimator of the observed
#					regressors (beta.eup)
#	  3. \$opt.d			= optimal dimension according to the
#					dimensionalty criterion
#	  4. \$nbr.iterations 	= number of iteration
###########################################################################
Eup.default <- function(formula,
additive.effects = c("none", "individual", "time", "twoways"),
dim.criterion	 = c("PC1", "PC2", "PC3", "BIC3", "IC1", "IC2" , "IC3", "IPC1", "IPC2", "IPC3"),
d.max            = NULL,
sig2.hat         = NULL,
factor.dim       = NULL,
double.iteration = TRUE,
start.beta       = NULL,
max.iteration    = 500,
convergence      = 1e-6,
...){

### substruct data from fomrmula and perfome a transformation according

# check fomula

if(!class(formula)=="formula"){
stop("\n Argument >>formula<< needs a formula-object like
y~x1+... where the elements are matrices")
}

# names

names  <- names(model.frame(formula))

sig2.hat.dim <- sig2.hat

## Transform the response variable as well as the 'P' regressors and give them in a list where ech
## componente contains also a list with:
#       1- "Tr" Name of the transformation
#       2- "I" Logical variable if ther is intercept or no
#       3- "ODM" Original Data matrix
#       4- "TDM" Transformed Data in a matrix
#       5- "TDV" Transformed Data in a NT x 1 Vector
#       6- "TRm" Sublist with
#           a- "OVc" Overall Constant
#           b- "InC" time constant individual effects
#           c- "TiVC" additive time varying effects

PF.obj	<- FUN.Pformula(formula = formula
nc 		<- ncol(PF.obj[[1]]\$ODM)
nr 		<- nrow(PF.obj[[1]]\$ODM)
P  		<- length(PF.obj) - 1
intercept <- PF.obj[[1]]\$I

# prepare for FUN.default (give the data in a large matrix: (y, x1, ...xP) where y, xp are  NT x1 Vectors  )

dat.dim 	  <- c(nr, nc, P)
dat.matrix	  <- sapply(1:(P+1), function(i) PF.obj[[i]]\$TDV)

dim.criterion <- match.arg(dim.criterion)

# Estimation results
tr.model.est <- FUN.Eup(dat.matrix		= dat.matrix
, dat.dim		= dat.dim
, double.iteration= double.iteration
, start.beta	= start.beta
, dim.criterion 	= dim.criterion
, factor.dim	= factor.dim
, d.max		= d.max
, max.iteration   = max.iteration
, convergence     = convergence
, sig2.hat = sig2.hat)

# Eup beta and Nbr.iteration
converges         <- tr.model.est\$converges
Nbr.iteration	<- tr.model.est\$nbr.iterations
beta.Eup		<- as.matrix(tr.model.est\$beta)
colnames(beta.Eup) <- ""
rownames(beta.Eup) <- names[-1]

# calculate the constant
OvConst <- sapply(1:(P+1),  function(i) PF.obj[[i]]\$TRm\$OVc)
ConsCoef <- OvConst[1] - OvConst[-1]%*%beta.Eup

# calculate the additive indivudal effects

ind.means <- sapply(1:(P+1),  function(i) PF.obj[[i]]\$TRm\$InC)
Ind.Eff <- ind.means[,1, drop = FALSE] - ind.means[,-1, drop = FALSE]%*%beta.Eup - c(ConsCoef)

# calculate the additive time effecs

tim.means <- sapply(1:(P+1),  function(i) PF.obj[[i]]\$TRm\$TiVC)
Tim.Eff <- tim.means[,1, drop = FALSE] - tim.means[,-1, drop = FALSE]%*%beta.Eup - c(ConsCoef)

# factor dimension
used.dim <- tr.model.est\$opt.d
proposed.dim <- ifelse(is.null(factor.dim), "Not Indicated"
, factor.dim)
if(is.null(factor.dim)){
optimal.dim <- tr.model.est\$opt.d
}
else{
optimal.dim <- EstDim(tr.model.est\$PCA
, dim.criterion = dim.criterion
, d.max = d.max, sig2.hat= sig2.hat)[,2]
}

# factor structure and resuduals

restrict.fs.a.resid	<- restrict.pca(tr.model.est\$PCA
, restrict.mode=restrict.mode)
fs.and.resid		<- restrict.fs.a.resid\$orig.values
factors			<- restrict.fs.a.resid\$factors
unob.fact.stru		<- restrict.fs.a.resid\$fitted.values
residuals			<- fs.and.resid - unob.fact.stru

# fitted values

orig.Y <- PF.obj[[1]]\$ODM
fitted.values <- orig.Y - residuals
degrees.of.freedom <- (nr*nc - (nr+nc)*used.dim - P -
intercept -
#sig2.hat <- sum(diag(crossprod(residuals)))/degrees.of.freedom
sig2.hat <- sum(diag(var(residuals)))*(nr-1)/degrees.of.freedom
## Results

final.result <- list(
dat.matrix = dat.matrix
, formula = formula
, dat.dim = dat.dim
, slope.para = beta.Eup
, names = names
, is.intercept = intercept
, Intercept = c(ConsCoef)
, unob.factors = factors
, unob.fact.stru = unob.fact.stru
, used.dim= used.dim
, proposed.dim= proposed.dim
, optimal.dim = optimal.dim
, factor.dim = factor.dim
, d.max = d.max
, dim.criterion = dim.criterion
, OvMeans = OvConst
, ColMean = ind.means
, RowMean = tim.means
, max.iteration = max.iteration
, convergence = convergence
, converges = converges
, start.beta = start.beta
, Nbr.iteration= Nbr.iteration
, fitted.values = fitted.values
, orig.Y = orig.Y
, residuals = residuals
, sig2.hat.dim = sig2.hat.dim
, sig2.hat = sig2.hat
, degrees.of.freedom = degrees.of.freedom
, restrict.mode = restrict.mode
, call = match.call())
class(final.result) <- "Eup"
return(final.result)
}
```

## Try the phtt package in your browser

Any scripts or data that you put into this service are public.

phtt documentation built on May 31, 2017, 4:17 a.m.