diff --git a/.github/workflows/cran-check.yml b/.github/workflows/cran-check.yml
index 404cb3c..326bbfc 100644
--- a/.github/workflows/cran-check.yml
+++ b/.github/workflows/cran-check.yml
@@ -22,7 +22,6 @@ jobs:
         uses: r-lib/actions/setup-r@v2
         with:
           r-version: 'release'
-          rtools-version: '42'
           update-rtools: true
       - name: Install R Dependencies
         uses: r-lib/actions/setup-r-dependencies@v2
diff --git a/NAMESPACE b/NAMESPACE
index 7e0b4f6..0f21cf3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -10,8 +10,10 @@ S3method(plot,gtdrm)
 S3method(plot,gwrm)
 S3method(plot,gwrmultiscalem)
 S3method(plot,modelselcritl)
+S3method(predict,gtwrm)
 S3method(predict,gwrm)
 S3method(print,gtdrm)
+S3method(print,gtwrm)
 S3method(print,gwrm)
 S3method(print,gwrmultiscalem)
 S3method(residuals,gtdrm)
@@ -22,6 +24,7 @@ S3method(step,gtdrm)
 S3method(step,gwrm)
 export(gtdr)
 export(gtdr_config)
+export(gtwr)
 export(gwr_basic)
 export(gwr_multiscale)
 export(mgwr_config)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 6c1c586..aaa6437 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,10 +1,18 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+gtwr_fit <- function(x, y, coords, times, bw, adaptive, kernel, lambda, longlat, p, theta, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, optim_lambda, verbose) {
+    .Call(`_GWmodel3_gtwr_fit`, x, y, coords, times, bw, adaptive, kernel, lambda, longlat, p, theta, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, optim_lambda, verbose)
+}
+
 gtdr_fit <- function(x, y, coords, bw, adaptive, kernel, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, optim_threashold, optim_step, optim_max_iter, select_model, select_model_threshold, variable_names, verbose) {
     .Call(`_GWmodel3_gtdr_fit`, x, y, coords, bw, adaptive, kernel, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, optim_threashold, optim_step, optim_max_iter, select_model, select_model_threshold, variable_names, verbose)
 }
 
+gtwr_fit <- function(x, y, coords, times, bw, adaptive, kernel, lambda, longlat, p, theta, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, optim_lambda, verbose) {
+    .Call(`_GWmodel3_gtwr_fit`, x, y, coords, times, bw, adaptive, kernel, lambda, longlat, p, theta, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, optim_lambda, verbose)
+}
+
 gwr_basic_fit <- function(x, y, coords, bw, adaptive, kernel, longlat, p, theta, optim_bw_lower, optim_bw_upper, hatmatrix, intercept, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, select_model, select_model_criterion, select_model_threshold, variable_names, verbose) {
     .Call(`_GWmodel3_gwr_basic_fit`, x, y, coords, bw, adaptive, kernel, longlat, p, theta, optim_bw_lower, optim_bw_upper, hatmatrix, intercept, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, select_model, select_model_criterion, select_model_threshold, variable_names, verbose)
 }
diff --git a/R/gtwr.R b/R/gtwr.R
new file mode 100644
index 0000000..7b5d10b
--- /dev/null
+++ b/R/gtwr.R
@@ -0,0 +1,292 @@
+#' Geographically and Temporally Weighted Regression
+#'
+#' @param formula Formula for regression.
+#' @param data A `sf` objects.
+#' @param timestamps A vector timestamps for all samples.
+#'  Either a numerical vector, or a character vector to be parsed according to `time_format`.
+#' @param bw Either a value to set the size of bandwidth,
+#'  or one of the following characters to set the criterion for
+#'  bandwidth auto-optimization process.
+#'  - `AIC`
+#'  - `CV`
+#'  Note that if `NA` or other non-numeric value is setted,
+#'  this parameter will be reset to `Inf`.
+#' @param lambda Either a value between 0 and 1 as the weight of temporal distance,
+#'  or a `NA` value to enable auto-selection.
+#' @param adaptive Whether the bandwidth value is adaptive or not.
+#' @param kernel Kernel function used.
+#' @param longlat Whether the coordinates
+#' @param p Power of the Minkowski distance,
+#'  default to 2, i.e., Euclidean distance.
+#' @param theta Angle in radian to roate the coordinate system, default to 0.
+#' @param optim_bw_range Bounds on bandwidth optimization, a vector of two numeric elements.
+#'  Set to `NA_real_` to enable default values selected by the algorithm.
+#' @param hatmatrix If TRUE, great circle will be caculated.
+#' @param time_format The format used to parse timestamps if they are characters.
+#' @param parallel_method Parallel method.
+#' @param parallel_arg Parallel method argument.
+#' @param verbose Whether to print additional information.
+#'
+#' @return A `gtwrm` object.
+#' 
+#' @details
+#' ## Parallelization
+#' 
+#' Supported method(s):
+#' 
+#' - Multithreading (`omp`)
+#' 
+#' See the vignettes about parallelization to learn more about this topic.
+#'
+#' @examples
+#' data(LondonHP)
+#' LondonHP$time <- as.integer(round(runif(nrow(LondonHP), 1, 365)))
+#'
+#' # Basic usage
+#' m <- gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, "time", 64, 0.05, TRUE)
+#' m
+#' head(coef(m))
+#' head(fitted(m))
+#' head(residuals(m))
+#'
+#' # Bandwidth Optimization
+#' gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, NA, 0.05, TRUE)
+#' 
+#' # Lambda Optimization
+#' gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, 64, NA, TRUE)
+#' 
+#' # Bandwidth and Lambda optimization
+#' gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, NA, NA, TRUE)
+#'
+#' @importFrom stats na.action model.frame model.extract model.matrix terms
+#' @export
+gtwr <- function(
+    formula,
+    data,
+    timestamps,
+    bw = NA,
+    lambda = NA,
+    adaptive = FALSE,
+    kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"),
+    longlat = FALSE,
+    p = 2.0,
+    theta = 0.0,
+    optim_bw_range = c(0, Inf),
+    hatmatrix = TRUE,
+    time_format = NA,
+    parallel_method = c("no", "omp"),
+    parallel_arg = c(0),
+    verbose = FALSE
+) {
+    ### Check args
+    kernel = match.arg(kernel)
+    parallel_method = match.arg(parallel_method)
+    attr(data, "na.action") <- getOption("na.action")
+
+    ### Extract coords
+    data <- do.call(na.action(data), args = list(data))
+    coords <- as.matrix(sf::st_coordinates(sf::st_centroid(data)))
+    if (is.null(coords) || nrow(coords) != nrow(data))
+        stop("Missing coordinates.")
+
+    ### Extract variables
+    mc <- match.call(expand.dots = FALSE)
+    mf <- model.frame(formula = formula(formula), data = sf::st_drop_geometry(data))
+    mt <- attr(mf, "terms")
+    y <- model.extract(mf, "response")
+    x <- model.matrix(mt, mf)
+    dep_var <- as.character(attr(terms(mf), "variables")[[2]])
+    has_intercept <- attr(terms(mf), "intercept") == 1
+    indep_vars <- colnames(x)
+    indep_vars[which(indep_vars == "(Intercept)")] <- "Intercept"
+    colnames(x) <- indep_vars
+    if (has_intercept && indep_vars[1] != "Intercept") {
+        stop("Please put Intercept to the first column.")
+    }
+
+    if (missing(timestamps)) {
+        stop("Timestamps are required.")
+    } else if (is.character(timestamps)) {
+        if (length(timestamps) == 1) {
+            if (timestamps %in% names(data)) {
+                timestamps <- data[[timestamps]]
+            } else {
+                stop("The specified `timestamps` must be a (only one) valid column name in `data`.")
+            }
+        } else if (length(timestamps) == nrow(coords)) {
+            timestamps <- as.numeric(as.Date(timestamps, format = time_format))
+        } else {
+            stop("If `timestamps` is characters, it must be either a column name in `data`, or date-time character vector for all samples.")
+        }
+    } else if (is.numeric(timestamps) || is.integer(timestamps)) {
+        if (length(timestamps) != nrow(coords)) {
+            stop("The length of `timestamps` must be equal to the number of samples.")
+        }
+    }
+
+    ### Check whether bandwidth is valid.
+    if (missing(bw)) {
+        optim_bw <- TRUE
+        optim_bw_criterion <- "AIC"
+        bw <- Inf
+    } else if (is.numeric(bw) || is.integer(bw)) {
+        optim_bw <- FALSE
+        optim_bw_criterion <- "AIC"
+    } else {
+        optim_bw <- TRUE
+        optim_bw_criterion <-
+            ifelse(is.character(bw), match.arg(bw, c("CV", "AIC")), "AIC")
+        bw <- Inf
+    }
+
+    if (missing(lambda)) {
+        optim_lambda <- TRUE
+        lambda <- 0.05
+    } else if (is.numeric(lambda) && lambda >= 0 && lambda <= 1) {
+        optim_lambda <- FALSE
+    } else {
+        optim_lambda <- TRUE
+        lambda <- 0.05
+    }
+
+    ### Call solver
+    c_result <- tryCatch(gtwr_fit(
+        x, y, coords, timestamps, bw, adaptive, enum(kernel), lambda, longlat, p, theta,
+        has_intercept, hatmatrix,
+        enum_list(parallel_method, parallel_types), parallel_arg,
+        optim_bw, enum(optim_bw_criterion, c("AIC", "CV")),
+        optim_lambda, as.integer(verbose)
+    ), error = function (e) {
+        stop("Error:", conditionMessage(e))
+    })
+    if (optim_bw)
+        bw <- c_result$bandwidth
+    betas <- c_result$betas
+    betas_se <- c_result$betasSE
+    shat_trace <- c_result$sTrace
+    fitted <- c_result$fitted
+    diagnostic <- c_result$diagnostic
+    resi <- y - fitted
+    n_dp <- nrow(coords)
+    rss_gw <- sum(resi * resi)
+    sigma <- rss_gw / (n_dp - 2 * shat_trace[1] + shat_trace[2])
+    betas_se <- sqrt(sigma * betas_se)
+    betas_tv <- betas / betas_se
+
+    ### Create result Layer
+    colnames(betas) <- indep_vars
+    colnames(betas_se) <- paste(indep_vars, "SE", sep = ".")
+    colnames(betas_tv) <- paste(indep_vars, "TV", sep = ".")
+    sdf_data <- as.data.frame(cbind(
+        betas,
+        "yhat" = fitted,
+        "residual" = resi,
+        betas_se,
+        betas_tv
+    ))
+    sdf_data$geometry <- sf::st_geometry(data)
+    sdf <- sf::st_sf(sdf_data)
+
+    ### Return result
+    gtwrm <- list(
+        SDF = sdf,
+        diagnostic = diagnostic,
+        args = list(
+            x = x,
+            y = y,
+            coords = coords,
+            timestamps = timestamps,
+            bw = bw,
+            adaptive = adaptive,
+            kernel = kernel,
+            longlat = longlat,
+            p = p,
+            theta = theta,
+            hatmatrix = hatmatrix,
+            has_intercept = has_intercept,
+            parallel_method = parallel_method,
+            parallel_arg = parallel_arg,
+            optim_bw = optim_bw,
+            optim_bw_criterion = optim_bw_criterion,
+            optim_lambda = optim_lambda,
+            verbose = verbose
+        ),
+        call = mc,
+        indep_vars = indep_vars,
+        dep_var = dep_var
+    )
+    class(gtwrm) <- c("gtwrm", "gwrm")
+    gtwrm
+}
+
+#' @describeIn print.gwrm
+#'
+#' @param x Object returned by GW modelling methods.
+#' @param decimal_fmt The format string passing to [base::sprintf()].
+#' @inheritDotParams print_table_md
+#' 
+#' @method print gtwrm
+#' @importFrom stats coef fivenum
+#' @rdname print
+#' @export
+print.gtwrm <- function(x, decimal_fmt = "%.3f", ...) {
+    if (!inherits(x, "gtwrm")) {
+        stop("It's not a gtwrm object.")
+    }
+
+    ### Basic Information
+    cat("Geographically and Temporally Weighted Regression Model", fill = T)
+    cat("=======================================================", fill = T)
+    cat("  Formula:", deparse(x$call$formula), fill = T)
+    cat("     Data:", deparse(x$call$data), fill = T)
+    cat("   Kernel:", x$args$kernel, fill = T)
+    cat("Bandwidth:", x$args$bw,
+        ifelse(x$args$adaptive, "(Nearest Neighbours)", "(Meters)"),
+        ifelse(x$args$optim_bw, paste0(
+            "(Optimized accroding to ",
+            x$args$optim_bw_criterion,
+            ")"
+        ), ""), fill = T)
+    cat("\n", fill = T)
+
+    cat("Summary of Coefficient Estimates", fill = T)
+    cat("--------------------------------", fill = T)
+    betas <- coef(x)
+    beta_fivenum <- t(apply(betas, 2, fivenum))
+    colnames(beta_fivenum) <- c("Min.", "1st Qu.", "Median", "3rd Qu.", "Max.")
+    rownames(beta_fivenum) <- colnames(betas)
+    beta_str <- rbind(
+        c("Coefficient", colnames(beta_fivenum)),
+        cbind(rownames(beta_fivenum), matrix2char(beta_fivenum, decimal_fmt))
+    )
+    print_table_md(beta_str, ...)
+    cat("\n", fill = T)
+
+    cat("Diagnostic Information", fill = T)
+    cat("----------------------", fill = T)
+    cat("  RSS:", x$diagnostic$RSS, fill = T)
+    cat("  ENP:", x$diagnostic$ENP, fill = T)
+    cat("  EDF:", x$diagnostic$EDF, fill = T)
+    cat("   R2:", x$diagnostic$RSquare, fill = T)
+    cat("R2adj:", x$diagnostic$RSquareAdjust, fill = T)
+    cat("  AIC:", x$diagnostic$AIC, fill = T)
+    cat(" AICc:", x$diagnostic$AICc, fill = T)
+    cat("\n", fill = T)
+}
+
+#' @describeIn gtwr [Not implemented] Predict on new locations.
+#'
+#' @param object A "gtwrm" object.
+#' @param regression_points Data of new locations.
+#' @param \dots Additional arguments.
+#' @param verbose Whether to print additional message.
+#' 
+#' @method predict gtwrm
+#' @noRd
+#' @export
+predict.gtwrm <- function(object, regression_points, verbose = FALSE, ...) {
+    if (!inherits(object, "gtwrm")) {
+        stop("It's not a gtwrm object.")
+    }
+    stop("GTWR is not predictable now.")
+}
diff --git a/R/gwr_basic.R b/R/gwr_basic.R
index b5a48c2..53f1864 100644
--- a/R/gwr_basic.R
+++ b/R/gwr_basic.R
@@ -311,9 +311,9 @@ step.gwrm <- function(
     object
 }
 
-#' Print description of a `gwrm` object
+#' Print description of a GW model object
 #'
-#' @param x An `hgwrm` object returned by [gwr_basic()].
+#' @param x Object returned by GW modelling methods.
 #' @param decimal_fmt The format string passing to [base::sprintf()].
 #' @inheritDotParams print_table_md
 #' 
diff --git a/R/gwr_multiscale.R b/R/gwr_multiscale.R
index 9b6b5e5..77e8b66 100644
--- a/R/gwr_multiscale.R
+++ b/R/gwr_multiscale.R
@@ -268,9 +268,9 @@ gwr_multiscale <- function(
     gwrmultiscalem
 }
 
-#' Print description of a `gwrmultiscalem` object
+#' @describeIn print.gwrm
 #'
-#' @param x An `gwrmultiscalem` object returned by [gwr_multiscale()].
+#' @param x Object returned by GW modelling methods.
 #' @param decimal_fmt The format string passing to [base::sprintf()].
 #' @inheritDotParams print_table_md
 #'
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 7a9b8d5..3db4e53 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -24,9 +24,10 @@ reference:
   - contents:
       - GWmodel3-package
   - title: GW regression models
-  - subtitle: Basic GWR
+  - subtitle: Basic GWR and Extensions
   - contents:
     - starts_with("gwr_basic")
+    - starts_with("gtwr")
   - subtitle: Multiscale GWR
   - contents:
     - starts_with("gwr_multiscale")
diff --git a/man/gtwr.Rd b/man/gtwr.Rd
new file mode 100644
index 0000000..0a51c08
--- /dev/null
+++ b/man/gtwr.Rd
@@ -0,0 +1,108 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gtwr.R
+\name{gtwr}
+\alias{gtwr}
+\title{Geographically and Temporally Weighted Regression}
+\usage{
+gtwr(
+  formula,
+  data,
+  timestamps,
+  bw = NA,
+  lambda = NA,
+  adaptive = FALSE,
+  kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"),
+  longlat = FALSE,
+  p = 2,
+  theta = 0,
+  optim_bw_range = c(0, Inf),
+  hatmatrix = TRUE,
+  time_format = NA,
+  parallel_method = c("no", "omp"),
+  parallel_arg = c(0),
+  verbose = FALSE
+)
+}
+\arguments{
+\item{formula}{Formula for regression.}
+
+\item{data}{A \code{sf} objects.}
+
+\item{timestamps}{A vector timestamps for all samples.
+Either a numerical vector, or a character vector to be parsed according to \code{time_format}.}
+
+\item{bw}{Either a value to set the size of bandwidth,
+or one of the following characters to set the criterion for
+bandwidth auto-optimization process.
+\itemize{
+\item \code{AIC}
+\item \code{CV}
+Note that if \code{NA} or other non-numeric value is setted,
+this parameter will be reset to \code{Inf}.
+}}
+
+\item{lambda}{Either a value between 0 and 1 as the weight of temporal distance,
+or a \code{NA} value to enable auto-selection.}
+
+\item{adaptive}{Whether the bandwidth value is adaptive or not.}
+
+\item{kernel}{Kernel function used.}
+
+\item{longlat}{Whether the coordinates}
+
+\item{p}{Power of the Minkowski distance,
+default to 2, i.e., Euclidean distance.}
+
+\item{theta}{Angle in radian to roate the coordinate system, default to 0.}
+
+\item{optim_bw_range}{Bounds on bandwidth optimization, a vector of two numeric elements.
+Set to \code{NA_real_} to enable default values selected by the algorithm.}
+
+\item{hatmatrix}{If TRUE, great circle will be caculated.}
+
+\item{time_format}{The format used to parse timestamps if they are characters.}
+
+\item{parallel_method}{Parallel method.}
+
+\item{parallel_arg}{Parallel method argument.}
+
+\item{verbose}{Whether to print additional information.}
+}
+\value{
+A \code{gtwrm} object.
+}
+\description{
+Geographically and Temporally Weighted Regression
+}
+\details{
+\subsection{Parallelization}{
+
+Supported method(s):
+\itemize{
+\item Multithreading (\code{omp})
+}
+
+See the vignettes about parallelization to learn more about this topic.
+}
+}
+\examples{
+data(LondonHP)
+LondonHP$time <- as.integer(round(runif(nrow(LondonHP), 1, 365)))
+
+# Basic usage
+m <- gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, "time", 64, 0.05, TRUE)
+m
+head(coef(m))
+head(fitted(m))
+head(residuals(m))
+
+# Bandwidth Optimization
+gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, NA, 0.05, TRUE)
+
+# Lambda Optimization
+gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, 64, NA, TRUE)
+
+# Bandwidth and Lambda optimization
+gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, NA, NA, TRUE)
+
+}
diff --git a/man/print.Rd b/man/print.Rd
index e7e67ee..a6233f8 100644
--- a/man/print.Rd
+++ b/man/print.Rd
@@ -1,24 +1,28 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/gtdr.R, R/gwr_basic.R, R/gwr_multiscale.R
+% Please edit documentation in R/gtdr.R, R/gtwr.R, R/gwr_basic.R,
+%   R/gwr_multiscale.R
 \name{print.gtdrm}
 \alias{print.gtdrm}
+\alias{print.gtwrm}
 \alias{print.gwrm}
 \alias{print.gwrmultiscalem}
 \title{Print description of a \code{gtdrm} object}
 \usage{
 \method{print}{gtdrm}(x, decimal_fmt = "\%.3f", ...)
 
+\method{print}{gtwrm}(x, decimal_fmt = "\%.3f", ...)
+
 \method{print}{gwrm}(x, decimal_fmt = "\%.3f", ...)
 
 \method{print}{gwrmultiscalem}(x, decimal_fmt = "\%.3f", ...)
 }
 \arguments{
-\item{x}{An \code{gwrmultiscalem} object returned by \code{\link[=gwr_multiscale]{gwr_multiscale()}}.}
+\item{x}{Object returned by GW modelling methods.}
 
 \item{decimal_fmt}{The format string passing to \code{\link[base:sprintf]{base::sprintf()}}.}
 
 \item{...}{
-  Arguments passed on to \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}
+  Arguments passed on to \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}
   \describe{
     \item{\code{col.sep}}{Column seperator. Default to \code{""}.}
     \item{\code{header.sep}}{Header seperator. Default to \code{"-"}.}
@@ -33,7 +37,5 @@ Possible values are \code{"plain"}, \code{"md"} or \code{"latex"}. Default to \c
 \description{
 Print description of a \code{gtdrm} object
 
-Print description of a \code{gwrm} object
-
-Print description of a \code{gwrmultiscalem} object
+Print description of a GW model object
 }
diff --git a/src/Makevars.in b/src/Makevars.in
index f5897d0..9806eed 100644
--- a/src/Makevars.in
+++ b/src/Makevars.in
@@ -34,6 +34,7 @@ endif
 OBJECTS_LIBGWMODEL = \
 	libgwmodel/src/gwmodelpp/spatialweight/BandwidthWeight.o \
 	libgwmodel/src/gwmodelpp/spatialweight/CRSDistance.o \
+	libgwmodel/src/gwmodelpp/spatialweight/CRSSTDistance.o \
 	libgwmodel/src/gwmodelpp/spatialweight/DMatDistance.o \
 	libgwmodel/src/gwmodelpp/spatialweight/Distance.o \
 	libgwmodel/src/gwmodelpp/spatialweight/MinkwoskiDistance.o \
@@ -41,6 +42,8 @@ OBJECTS_LIBGWMODEL = \
 	libgwmodel/src/gwmodelpp/spatialweight/SpatialWeight.o \
 	libgwmodel/src/gwmodelpp/spatialweight/Weight.o \
 	libgwmodel/src/gwmodelpp/BandwidthSelector.o \
+	libgwmodel/src/gwmodelpp/GTWR.o \
+	libgwmodel/src/gwmodelpp/GWPCA.o \
 	libgwmodel/src/gwmodelpp/GTDR.o \
 	libgwmodel/src/gwmodelpp/GWRBase.o \
 	libgwmodel/src/gwmodelpp/GWRBasic.o \
@@ -60,6 +63,7 @@ OBJECTS_GWMODEL = \
 	utils.o \
 	gwr_basic.o \
 	gwr_multiscale.o \
+	gtwr.o \
 	gtdr.o \
 	RcppExports.o
 
diff --git a/src/Makevars.win b/src/Makevars.win
index 8ac438f..e4f42d4 100644
--- a/src/Makevars.win
+++ b/src/Makevars.win
@@ -36,6 +36,7 @@ PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(GSL_L
 OBJECTS_LIBGWMODEL = \
 	libgwmodel/src/gwmodelpp/spatialweight/BandwidthWeight.o \
 	libgwmodel/src/gwmodelpp/spatialweight/CRSDistance.o \
+	libgwmodel/src/gwmodelpp/spatialweight/CRSSTDistance.o \
 	libgwmodel/src/gwmodelpp/spatialweight/DMatDistance.o \
 	libgwmodel/src/gwmodelpp/spatialweight/Distance.o \
 	libgwmodel/src/gwmodelpp/spatialweight/MinkwoskiDistance.o \
@@ -43,6 +44,7 @@ OBJECTS_LIBGWMODEL = \
 	libgwmodel/src/gwmodelpp/spatialweight/SpatialWeight.o \
 	libgwmodel/src/gwmodelpp/spatialweight/Weight.o \
 	libgwmodel/src/gwmodelpp/BandwidthSelector.o \
+	libgwmodel/src/gwmodelpp/GTWR.o \
 	libgwmodel/src/gwmodelpp/GTDR.o \
 	libgwmodel/src/gwmodelpp/GWRBase.o \
 	libgwmodel/src/gwmodelpp/GWRBasic.o \
@@ -63,6 +65,7 @@ OBJECTS_GWMODEL = \
 	gwr_basic.o \
 	gwr_multiscale.o \
 	gtdr.o \
+	gtwr.o \
 	RcppExports.o
 
 
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index d0cedbc..2fb726d 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -40,6 +40,35 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// gtwr_fit
+List gtwr_fit(const arma::mat& x, const arma::vec& y, const arma::mat& coords, const arma::vec& times, double bw, bool adaptive, int kernel, double lambda, bool longlat, double p, double theta, bool intercept, bool hatmatrix, size_t parallel_type, const IntegerVector& parallel_arg, bool optim_bw, size_t optim_bw_criterion, bool optim_lambda, int verbose);
+RcppExport SEXP _GWmodel3_gtwr_fit(SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP timesSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP lambdaSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP interceptSEXP, SEXP hatmatrixSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP, SEXP optim_bwSEXP, SEXP optim_bw_criterionSEXP, SEXP optim_lambdaSEXP, SEXP verboseSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type coords(coordsSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type times(timesSEXP);
+    Rcpp::traits::input_parameter< double >::type bw(bwSEXP);
+    Rcpp::traits::input_parameter< bool >::type adaptive(adaptiveSEXP);
+    Rcpp::traits::input_parameter< int >::type kernel(kernelSEXP);
+    Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
+    Rcpp::traits::input_parameter< bool >::type longlat(longlatSEXP);
+    Rcpp::traits::input_parameter< double >::type p(pSEXP);
+    Rcpp::traits::input_parameter< double >::type theta(thetaSEXP);
+    Rcpp::traits::input_parameter< bool >::type intercept(interceptSEXP);
+    Rcpp::traits::input_parameter< bool >::type hatmatrix(hatmatrixSEXP);
+    Rcpp::traits::input_parameter< size_t >::type parallel_type(parallel_typeSEXP);
+    Rcpp::traits::input_parameter< const IntegerVector& >::type parallel_arg(parallel_argSEXP);
+    Rcpp::traits::input_parameter< bool >::type optim_bw(optim_bwSEXP);
+    Rcpp::traits::input_parameter< size_t >::type optim_bw_criterion(optim_bw_criterionSEXP);
+    Rcpp::traits::input_parameter< bool >::type optim_lambda(optim_lambdaSEXP);
+    Rcpp::traits::input_parameter< int >::type verbose(verboseSEXP);
+    rcpp_result_gen = Rcpp::wrap(gtwr_fit(x, y, coords, times, bw, adaptive, kernel, lambda, longlat, p, theta, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw, optim_bw_criterion, optim_lambda, verbose));
+    return rcpp_result_gen;
+END_RCPP
+}
 // gwr_basic_fit
 List gwr_basic_fit(const arma::mat& x, const arma::vec& y, const arma::mat& coords, double bw, bool adaptive, size_t kernel, bool longlat, double p, double theta, double optim_bw_lower, double optim_bw_upper, bool hatmatrix, bool intercept, size_t parallel_type, const IntegerVector& parallel_arg, bool optim_bw, size_t optim_bw_criterion, bool select_model, size_t select_model_criterion, size_t select_model_threshold, const CharacterVector& variable_names, int verbose);
 RcppExport SEXP _GWmodel3_gwr_basic_fit(SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP optim_bw_lowerSEXP, SEXP optim_bw_upperSEXP, SEXP hatmatrixSEXP, SEXP interceptSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP, SEXP optim_bwSEXP, SEXP optim_bw_criterionSEXP, SEXP select_modelSEXP, SEXP select_model_criterionSEXP, SEXP select_model_thresholdSEXP, SEXP variable_namesSEXP, SEXP verboseSEXP) {
@@ -134,6 +163,7 @@ END_RCPP
 
 static const R_CallMethodDef CallEntries[] = {
     {"_GWmodel3_gtdr_fit", (DL_FUNC) &_GWmodel3_gtdr_fit, 19},
+    {"_GWmodel3_gtwr_fit", (DL_FUNC) &_GWmodel3_gtwr_fit, 19},
     {"_GWmodel3_gwr_basic_fit", (DL_FUNC) &_GWmodel3_gwr_basic_fit, 22},
     {"_GWmodel3_gwr_basic_predict", (DL_FUNC) &_GWmodel3_gwr_basic_predict, 14},
     {"_GWmodel3_gwr_multiscale_fit", (DL_FUNC) &_GWmodel3_gwr_multiscale_fit, 25},
diff --git a/src/gtwr.cpp b/src/gtwr.cpp
new file mode 100644
index 0000000..a60b087
--- /dev/null
+++ b/src/gtwr.cpp
@@ -0,0 +1,124 @@
+// [[Rcpp::depends(RcppArmadillo)]]
+#include <RcppArmadillo.h>
+#include "utils.h"
+#include "gwmodelpp/GTWR.h"
+
+using namespace std;
+using namespace Rcpp;
+using namespace arma;
+using namespace gwm;
+
+// [[Rcpp::export]]
+List gtwr_fit(
+    const arma::mat& x,
+    const arma::vec& y,
+    const arma::mat& coords,
+    const arma::vec& times,
+    double bw,
+    bool adaptive,
+    int kernel,
+    double lambda,
+    bool longlat,
+    double p,
+    double theta,
+    bool intercept,
+    bool hatmatrix,
+    size_t parallel_type,
+    const IntegerVector& parallel_arg,
+    bool optim_bw,
+    size_t optim_bw_criterion,
+    bool optim_lambda,
+    int verbose
+) {
+    // Convert data types
+    std::vector<int> vpar_args = as< std::vector<int> >(Rcpp::IntegerVector(parallel_arg));
+
+    // Make Spatial Weight
+    BandwidthWeight bandwidth(bw, adaptive, BandwidthWeight::KernelFunctionType(kernel));
+    CRSDistance *sdist = nullptr;
+    if (longlat)
+    {
+        sdist = new CRSDistance(true);
+    }
+    else
+    {
+        if (p == 2.0 && theta == 0.0)
+        {
+            sdist = new CRSDistance(false);
+        }
+        else
+        {
+            sdist = new MinkwoskiDistance(p, theta);
+        }
+    }
+    OneDimDistance tdist;
+    CRSSTDistance stdist(sdist, &tdist, lambda);
+    SpatialWeight spatial(&bandwidth, &stdist);
+    
+    // Make Algorithm Object
+    GTWR algorithm;
+    algorithm.setDependentVariable(y);
+    algorithm.setIndependentVariables(x);
+    algorithm.setCoords(coords, times);
+    algorithm.setSpatialWeight(spatial);
+    algorithm.setHasHatMatrix(hatmatrix);
+    algorithm.setHasIntercept(intercept);
+
+    if (optim_bw) 
+    {
+        algorithm.setIsAutoselectBandwidth(true);
+        algorithm.setBandwidthSelectionCriterion(GTWR::BandwidthSelectionCriterionType(size_t(optim_bw_criterion)));
+    }
+
+    if (optim_lambda)
+    {
+        algorithm.setIsAutoselectLambda(true);
+    }
+
+    switch (ParallelType(size_t(parallel_type)))
+    {
+    case ParallelType::SerialOnly:
+        algorithm.setParallelType(ParallelType::SerialOnly);
+        break;
+#ifdef _OPENMP
+    case ParallelType::OpenMP:
+        algorithm.setParallelType(ParallelType::OpenMP);
+        algorithm.setOmpThreadNum(vpar_args[0]);
+        break;
+#endif
+    default:
+        algorithm.setParallelType(ParallelType::SerialOnly);
+        break;
+    }
+    try
+    {
+        algorithm.fit();
+    }
+    catch(const std::exception& e)
+    {
+        stop(e.what());
+    }
+        
+    // Return Results
+    mat betas = algorithm.betas();
+    vec fitted = sum(x % betas, 1);
+    List result_list = List::create(
+        Named("betas") = betas,
+        Named("betasSE") = algorithm.betasSE(),
+        Named("sTrace") = algorithm.sHat(),
+        Named("sHat") = algorithm.s(),
+        Named("fitted") = fitted,
+        Named("diagnostic") = mywrap(algorithm.diagnostic())
+    );
+    const SpatialWeight& spatialWeights = algorithm.spatialWeight();
+    if (optim_bw)
+    {
+        result_list["bw_value"] = wrap(spatialWeights.weight<BandwidthWeight>()->bandwidth());
+    }
+    // if (optim_lambda)
+    // {
+    //     result_list["lambda_value"] = wrap(spatialWeights.distance<CRSSTDistance>()->mLambda);
+    // }
+
+    return result_list;
+}
diff --git a/tests/testthat/test-gtwr.R b/tests/testthat/test-gtwr.R
new file mode 100644
index 0000000..6a6e2fa
--- /dev/null
+++ b/tests/testthat/test-gtwr.R
@@ -0,0 +1,35 @@
+data(LondonHP)
+LondonHP$time <- as.integer(round(runif(nrow(LondonHP), 1, 365)))
+m <- NULL
+
+test_that("GTWR: set timestamps", {
+  m <<- expect_no_error({
+    gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, "time", 64, 0.05, TRUE)
+  })
+  expect_no_error({
+    gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, 64, 0.05, TRUE)
+  })
+})
+
+test_that("GTWR: optim lambda and bandwidth", {
+  expect_no_error({
+    gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, NA, 0.05, TRUE)
+  })
+  expect_no_error({
+    gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, 64, NA, TRUE)
+  })
+  expect_no_error({
+    gtwr(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, LondonHP$time, NA, NA, TRUE)
+  })
+})
+
+test_that("GTWR: helper functions", {
+  expect_no_error({
+    coef(m)
+    fitted(m)
+    residuals(m)
+  })
+  expect_no_error({
+    plot(m)
+  })
+})