I would like to replicate the p-value of NeweyWest function in the sandwich package. Could you let me know if there is any documentation for the explanation?

Sample dataset

# Question: R: Replicating P-value of `NeweyWest` of `sandwich` package.

library(sandwich)
library(lmtest)

# Sample dataset with 10 predictor variables
set.seed(123) 

n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
x4 <- rnorm(n)
x5 <- rnorm(n)
x6 <- rnorm(n)
x7 <- rnorm(n)
x8 <- rnorm(n)
x9 <- rnorm(n)
x10 <- rnorm(n)

y <- 2 + 1.5*x1 - 0.5*x2 + 0.3*x3 + 2*x4 + x5 - x6 + 0.8*x7 - 1.2*x8 + 0.5*x9 + 
  1.1*x10 + rnorm(n)

d <- data.frame(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, y)

eq <- y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10

Model

l <- lm(eq, d) 

# Create Newey-West covariance matrix estimates
nw_vcov <- NeweyWest(l, lag=1, prewhite=FALSE, adjust=TRUE)

# compute the square root of the diagonal elements in vcov
nw_se <- sqrt(diag(nw_vcov))

# Apply estimates to linear model
l_nw_se <- coeftest(l, vcov.=nw_vcov) 

# Extract the results
s <- as.data.frame(l_nw_se[, ])

# Question: how do I replicate the p-value? 
# Could you let me know if there is any documentation for the explanation?

P values in general are calculated using the Student t distribution pt() with t statistics and the degrees of freedom. For the usual two-tailed test, which is also applied in lmtest::coeftest, we can do:

> (p_values <- 2*pt(q=-abs(s$`t value`), df=l$df.residual))
 [1] 9.286402e-38 9.345499e-23 1.359084e-05 9.134816e-02 2.707676e-40
 [6] 2.393224e-15 7.685513e-12 5.457149e-13 1.957496e-15 6.150043e-04
[11] 5.182621e-21
> stopifnot(identical(p_values, s$`Pr(>|t|)`))

From sandwich::NeweyWest you are actually just creating the vcov where the standard errors are calculated from.