re-submission

cran
HujieBai · Apr 2, 2024 · ecd4f81 · ecd4f81
1 parent 62bc79b
commit ecd4f81
Show file tree

Hide file tree

Showing 19 changed files with 825 additions and 965 deletions.
diff --git a/.RData b/.RData
diff --git a/.Rhistory b/.Rhistory
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,16 +1,16 @@
 Package: DPTM
 Type: Package
-Title: Dynamic Panel Threshold Models with Fixed Effects Based on  Maximum Likelihood Estiamtion and MCMC
-Version: 1.2.3
+Title: Dynamic Panel Multiple Threshold Model with Fixed Effects
+Version: 1.3.5
 Date: 2024-03-28
 Authors@R: 
     person("Bai", "Hujie", , "[email protected]", role = c("aut", "cre", "cph"),
            comment = c(ORCID = "0009-0004-2060-4351"))
-Description: Methods and tools for dynamic panel data analysis. Allows estimation and model selection of linear model and threshold model.Multiple threshold model and model with time fixed effects are also allowed.
+Description:  Compute the fixed effects dynamic panel threshold model suggested by Ramírez-Rondán (2020) <doi:10.1080/07474938.2019.1624401>, and dynamic panel linear model suggested by Hsiao et al. (2002) <doi:10.1016/S0304-4076(01)00143-9>, where maximum likelihood type estimators are used. Multiple threshold estimation based on Markov Chain Monte Carlo (MCMC) is allowed, and model selection of linear model, threshold model and multiple threshold model is also allowed.
 License: GPL (>= 3)
 URL: https://github.com/HujieBai/DPTM
 Encoding: UTF-8
-Imports: Rcpp (>= 1.0.12),BayesianTools, purrr, MASS,stats,coda,parabar,utils
+Imports: Rcpp (>= 1.0.12),BayesianTools, purrr, MASS,stats,coda,parabar
 LinkingTo: Rcpp,RcppEigen
 RoxygenNote: 7.3.1
 Depends: 

diff --git a/NAMESPACE b/NAMESPACE
@@ -10,5 +10,4 @@ import(stats)
 importFrom(MASS,ginv)
 importFrom(coda,gelman.diag)
 importFrom(purrr,map_dbl)
-importFrom(utils,capture.output)
 useDynLib(DPTM)
diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,6 @@
+# DPTM 1.3.5
+* CRAN resubmission.
+
 # DPTM 1.2.3
 
 # DPTM 1.1.7

diff --git a/R/DPM.R b/R/DPM.R
@@ -1,6 +1,6 @@
 
-MLE <- function(y,x=NULL,delty0 =NULL,x1=NULL,cvs=NULL,ny=1,w=NULL,var_u = NULL,tt,nn,assumption = 1,
-                restart = FALSE,stages = 1){
+MLE <- function(y,x=NULL,delty0 =NULL,x1=NULL,cvs=NULL,ny=1,w=NULL,var_u = NULL,tt,nn,
+                restart = FALSE){
 
   if(is.null(delty0)){
     delty0 <- lag_transform(as.matrix(y),tt,nn,1,FALSE) -  lag_transform(as.matrix(y),tt,nn,1,TRUE)
@@ -42,12 +42,21 @@ MLE <- function(y,x=NULL,delty0 =NULL,x1=NULL,cvs=NULL,ny=1,w=NULL,var_u = NULL,
     stop("\n","The inital deltaX1 Matrix is singular!","\n",
          "Please change the parameter x1!","\n")
   }
-  delts <- qr.solve(x1int,as.vector(yint))
+
+
+
+
+  if(sum(x1int) == 0){
+    delts <- NULL
+  }else{
+    delts <- qr.solve(x1int,as.vector(yint))
+  }
 
   omega = diag(2,nrow = (tt-1),ncol = (tt-1))
   diag(omega[-1,-(tt-1)]) = -1
   diag(omega[-(tt-1),-1]) = -1
 
+
   deltxx <- cbind(deltxs,ininxc)
 
 
@@ -56,7 +65,7 @@ MLE <- function(y,x=NULL,delty0 =NULL,x1=NULL,cvs=NULL,ny=1,w=NULL,var_u = NULL,
   dyy <- as.matrix(deltxs[,1:ny])
   dxx <- deltxs[,(ny+1):ncol(deltxs)]
 
-  dx3 <-  lag_transform(dxx,tt-1,nn,2,FALSE)
+  dx3 <-  lag_transform(as.matrix(dxx),tt-1,nn,2,FALSE)
   dy1 <- cbind( lag_transform(dyy,tt-1,nn,2,FALSE),dx3)
   dy2 <- cbind( lag_transform( lag_transform(dyy,tt-1,nn,1,FALSE),tt-2,nn,1,TRUE),dx3)
 
@@ -81,70 +90,55 @@ MLE <- function(y,x=NULL,delty0 =NULL,x1=NULL,cvs=NULL,ny=1,w=NULL,var_u = NULL,
   varv1 <- (sum(yint^2)/nn)
 
 
-  if(assumption != 1){
-
-    if(!is.null(w)){
-      what <- w
+  if(!is.null(w)){
+    what <- w
+  }else{
+    what <- varv1/varu
+  }
+
+  if(what <= 1){
+    stop("\n","The initial w must > 1, thus please input a new w.","\n")
+  }
+
+  pars <- c(pars,what)
+  cd = length(pars)
+
+  result_in <- suppressWarnings(try(stats::nlm(three_two,pars,delty0=delty0,evs=deltxx,
+                                               omega=omega,cd=cd,tt=tt,nn=nn,iterlim = 500)
+                                    ,silent = TRUE))
+
+  if("try-error" %in% class(result_in)){
+    if(isTRUE(restart)){
+      pars[1:ny] = runif(ny,-1/ny,1/ny)
+      pars[cd-1] = runif(1,2,5)
+      pars[cd] = pars[cd-1]/varv1
+
+      result_in = suppressWarnings(try(stats::nlm(three_two,pars,delty0=delty0,evs=deltxx,
+                                                  omega=omega,cd=cd,tt=tt,nn=nn,iterlim = 500)
+                                       ,silent = TRUE))
     }else{
-      what <- varv1/varu
-    }
-
-    if(what <= 1){
-      stop("\n","When assumption == 2, initial w must > 1, thus please input a new w.","\n")
+      stop("\n","Encounter iteration failure!","\n",
+           "Please set restart = True or check other inputs!","\n")
     }
-
-    pars <- c(pars,what)
-    cd = length(pars)
-
-    result_in <- suppressWarnings(try(stats::nlm(three_two,pars,delty0=delty0,evs=deltxx,
-                                    omega=omega,cd=cd,tt=tt,nn=nn,iterlim = 500)
-                         ,silent = TRUE))
-
-    if("try-error" %in% class(result_in)){
-      if(isTRUE(restart)){
-        pars[1:ny] = runif(ny,-1/ny,1/ny)
-        pars[cd-1] = runif(1,2,5)
-        pars[cd] = pars[cd-1]/varv1
-
-        result_in = suppressWarnings(try(stats::nlm(three_two,pars,delty0=delty0,evs=deltxx,
-                                                    omega=omega,cd=cd,tt=tt,nn=nn,iterlim = 500)
-                                         ,silent = TRUE))
-      }else{
-        stop("\n","Encounter iteration failure!","\n",
-             "Please set restart = True or check other inputs!","\n")
-      }
-    }
-
-
-
-  }else{
-
-    losf <- ifelse(stages == 2,three_oneb,three_one)
-
-    cd = length(pars)
-
-    result_in = suppressWarnings(try(stats::nlm(losf,pars,delty0=delty0,evs=deltxx,
-                                                omega=omega,cd=cd,tt=tt,nn=nn,
-                                                ny=ny,varv1=varv1,iterlim = 500)
+  }
+
+  if(result_in$code == 1 & isTRUE(restart)){
+
+    pars[1:ny] = runif(ny,-1/ny,1/ny)
+    pars[cd-1] = runif(1,2,5)
+    pars[cd] = pars[cd-1]/varv1
+
+    result_in = suppressWarnings(try(stats::nlm(three_two,pars,delty0=delty0,evs=deltxx,
+                                                omega=omega,cd=cd,tt=tt,nn=nn,iterlim = 500)
                                      ,silent = TRUE))
-
-    if("try-error" %in% class(result_in)){
-      if(isTRUE(restart)){
-        pars[1:ny] = runif(ny,-1/ny,1/ny)
-        pars[cd] = varv1/2
-        result_in = suppressWarnings(try(stats::nlm(losf,pars,delty0=delty0,evs=deltxx,
-                                                    omega=omega,cd=cd,tt=tt,nn=nn,
-                                                    ny=ny,varv1=varv1,iterlim = 500)
-                                         ,silent = TRUE))
-      }else{
-        stop("\n","Encounter iteration failure!","\n",
-             "Please set restart = True or check other inputs!","\n")
-      }
-    }
-
-
-
+
   }
+
+  if("try-error" %in% class(result_in)){
+    stop("\n","Encounter iteration failure!","\n",
+         "Please set restart = True or check other inputs!","\n")
+  }
+
 
   pars <- result_in$estimate
   ssemin <- result_in$minimum
@@ -157,7 +151,7 @@ MLE <- function(y,x=NULL,delty0 =NULL,x1=NULL,cvs=NULL,ny=1,w=NULL,var_u = NULL,
 
   cms = matrix(0,nrow = evs,ncol = evs)
 
-  ww <- ifelse(assumption == 1,pars[cd-1],varv1/pars[cd])
+  ww <- varv1/pars[cd]
 
   xomega = omega
   xomega[1,1] = ww
@@ -211,28 +205,24 @@ MAP2 <- function (bayesianOutput, ...)
 #'@title The dynamic panel linear model with fixed effects
 #'@param y the dependent variable; vector type input.
 #'@param x the independent variable; matrix type input.
-#'@param y1 the lag dependent variable; vector type input; By default y1 is NULL,
+#'@param y1 the lag dependent variable; vector type input; By default, y1 is NULL,
 #'and then y1 will be computed by y automatically.
-#'@param time_trend the time trend; By default it is FALSE.
-#'@param time_fix_effects the time fixed effects; By default it is FALSE.
+#'@param time_trend the time trend; By default, it is FALSE.
+#'@param time_fix_effects the time fixed effects; By default, it is FALSE.
 #'@param x1 the initial values of independent variable; matrix type input.
-#'By default x1 is NULL, and thus x1 will be computed by x automatically.
+#'By default, x1 is NULL, and thus x1 will be computed by x automatically.
 #'@param tt the length of time period.
 #'@param nn the number of individuals.
-#'@param assumption the option of assumption; By default assumption is 1, and it can be 2;
-#'More details see Hsiao (2002).
-#'@param restart the option of iterations; By default restart is FALSE,
+#'@param restart the option of iterations; By default, restart is FALSE,
 #'if encounters iteration failure, please set restart as TRUE.
-#'@param Only_b the option of initial equation;By default Only_b is FALSE, and if Only_b is TRUE, initial delta y will be a constant C.;
-#'More details please see Hsiao (2002) and Ramirez-Rondan (2020).
-#'@param w the variance ratio; By default is NULL; It must be greater than 1, and
-#'only works when assumption is 2.
-#'@param var_u the option of variance of error term; By default is NULL; It must be
+#'@param w the variance ratio; By default, is NULL; It must be greater than 1.
+#'@param var_u the option of variance of error term; By default, is NULL; It must be
 #'greater than 0; When meet relevant ERROR, please change the var_u.
-#'@param delty0 the option of delta_y; By default delty0 is NULL; Pleas do not change delty0.
-#'@param Only_b the option of initial equation;By default Only_b is FALSE, and if Only_b is TRUE, initial delta y will be a constant C.
-#'More details please see Hsiao (2002) and Ramirez-Rondan (2020).
-#'@references Ramirez-Rondan, N. R. (2020). Maximum likelihood estimation
+#'@param delty0 the option of delta_y; By default, delty0 is NULL; Please do not change delty0.
+#'@param Only_b the option of initial equation;By default, Only_b is FALSE, and if Only_b is TRUE, initial delta y will be a constant C.
+#'Please see Hsiao (2002) and Ramírez-Rondán (2020) for more details.
+#'@param display the option of whether to print the messages of estimated results; By default, the display is TRUE.
+#'@references Ramírez-Rondán, N. R. (2020). Maximum likelihood estimation
 #' of dynamic panel threshold models. Econometric Reviews, 39(3), 260-276.
 #'@references Hsiao, C., Pesaran, M. H., & Tahmiscioglu, A. K. (2002).
 #' Maximum likelihood estimation of fixed effects dynamic panel data models covering short time periods. Journal of econometrics, 109(1), 107-150.
@@ -246,24 +236,34 @@ MAP2 <- function (bayesianOutput, ...)
 #'tt <- data$data_test_linear$tt
 #'nn <- data$data_test_linear$nn
 #'xx <- cbind(x,z)
-#'m1 <- DPML(y=y,x=xx,tt=tt,nn=nn,assumption = 1)
+#'m1 <- DPML(y=y,x=xx,tt=tt,nn=nn)
 #'m1$Coefs
 #'@describeIn DPML This is a dynamic panel linear model with fixed effects, which
 #'allows time trend term or time fixed effects.
-#'@returns A List of estimate results.
+#'@returns A list containing the following components:
+#'\item{ssemin}{  the negaive log-likelihood function value}
+#'\item{Coefs}{  parameter estimates containing t-values}
+#'\item{pars}{  iterated results for all parameters}
+#'\item{duit}{  the first-difference form of residuals}
+#'\item{dy0}{  the first-difference form of dependent variable}
+#'\item{xx}{  the independent variables and their initial values}
+#'\item{covariance_matrix}{  the covariance matrix}
+#'\item{Ses}{  the standard errors of coefs}
+#'\item{Zvalues}{  the values of the statistic}
+#'\item{ccd}{  the number of independent variables}
+#'\item{coefs}{  parameter estimates containing their initial valuess}
 #'@useDynLib DPTM
 #'@import Rcpp
 #'@import BayesianTools
 #'@import stats
 #'@import parabar
 #'@importFrom purrr map_dbl
-#'@importFrom utils capture.output
 #'@importFrom MASS ginv
 #'@importFrom coda gelman.diag
 #'@export
-DPML <- function(y,y1=NULL,x=NULL,w=NULL,var_u = NULL,tt,nn,assumption = 1,
+DPML <- function(y,y1=NULL,x=NULL,w=NULL,var_u = NULL,tt,nn,
                  time_trend =FALSE,time_fix_effects=FALSE,restart = FALSE,
-                 x1=NULL,delty0=NULL,Only_b = FALSE){
+                 x1=NULL,delty0=NULL,Only_b = FALSE,display = TRUE){
   ny=1
   time_shifts <- as.matrix(rep(1:tt,nn))
   time_effects <- kronecker(rep(1,nn),diag(tt))[,-c(1,2,3)]
@@ -296,7 +296,7 @@ DPML <- function(y,y1=NULL,x=NULL,w=NULL,var_u = NULL,tt,nn,assumption = 1,
   }
 
   fit_model <- MLE(y=y,x=cbind(y1,x),x1=x1,cvs=cvs0,ny=ny,w=w,var_u = var_u,tt=tt,nn=nn,
-             assumption = assumption,restart = restart,delty0=delty0)
+             restart = restart,delty0=delty0)
 
 
 
@@ -341,18 +341,21 @@ DPML <- function(y,y1=NULL,x=NULL,w=NULL,var_u = NULL,tt,nn,assumption = 1,
   colnames(jgs) <- c("Coefs","Significance","t-value")
 
   jgs <- jgs[1:fit_model$ccd,]
+
+  if(display == TRUE){
+    cat("\n","This is a Dynamic panel modedl with fixed effects.","\n")
+    cat("\n","---------------------------------------------------","\n")
+    cat("\n","Time Fixed Effects: ",time_fix_effects," !\n")
+    cat("\n","---------------------------------------------------","\n")
+    cat("\n","Time Shifts: ",time_trend," !\n")
+    cat("\n","---------------------------------------------------","\n")
+
+
+    cat("\n","The coefs are: ","\n")
+    print(jgs)
+
+  }
 
-  cat("\n","This is a Dynamic panel modedl with fixed effects.","\n",
-          "It follows Assumption ",assumption,", and Only_b =",Only_b,"!\n")
-  cat("\n","---------------------------------------------------","\n")
-  cat("\n","Time Fixed Effects: ",time_fix_effects," !\n")
-  cat("\n","---------------------------------------------------","\n")
-  cat("\n","Time Shifts: ",time_trend," !\n")
-  cat("\n","---------------------------------------------------","\n")
-
-
-  cat("\n","The coefs are: ","\n")
-  print(jgs)
 
 
 

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -9,14 +9,6 @@ lag_transform <- function(objs, t, n, lag, top) {
     .Call('_DPTM_lag_transform', PACKAGE = 'DPTM', objs, t, n, lag, top)
 }
 
-three_one <- function(pars1, delty0, evs, omega, cd, tt, nn, ny, varv1) {
-    .Call('_DPTM_three_one', PACKAGE = 'DPTM', pars1, delty0, evs, omega, cd, tt, nn, ny, varv1)
-}
-
-three_oneb <- function(pars1, delty0, evs, omega, cd, tt, nn, ny, varv1) {
-    .Call('_DPTM_three_oneb', PACKAGE = 'DPTM', pars1, delty0, evs, omega, cd, tt, nn, ny, varv1)
-}
-
 three_two <- function(pars1, delty0, evs, omega, cd, tt, nn) {
     .Call('_DPTM_three_two', PACKAGE = 'DPTM', pars1, delty0, evs, omega, cd, tt, nn)
 }