Determining derivatives from GAM smooth object

I have a quite simple time series data set consisting of annual averages of a singe variable ("AVERAGE"). I wish to investigate the rate of change (1st derivative) and acceleration (2nd derivative) and associated standard errors of the "trend" component of the time series. I have obtained the "trend" using the GAM and PREDICT functions of MGCV simply as follows:

A <- gam(AVERAGE ~ s(YEAR), data=DF, na.action=na.omit)
B <- predict(A, type="response", se.fit=TRUE)

I have determined derivatives through 2 separate methods, applying a high DoF cubic smooth spline and via first and second differences (lightly smoothed) and bootstrapping to approximate errors with both producing comparable results.

I note that the "gam.fit3" function facilitates determining upto 2nd order derivatives but is not called directly. I also note that using "predict.gam" with the type "lpmatrix" facilitates derivatives of smooths. I would like to use the "GAM" function directly to calculate the 1st and 2nd derivatives but am not sufficiently skilled to calculate or extract these derivatives. I tried to reconfigure Wood's example at the end of the "Predict.gam" help page for one variable but with no success. Any help to get me headed in the right direction would be terrific. Thanks Phil.

The example from predict.gam uses finite differences to approximate the derivatives of the smoothed terms

Here is an example to do this for a single predictor model. This is more straightforward that the example from the help.

A <- gam(AVERAGE ~ s(YEAR), data=DF, na.action=na.omit)
# new data for prediction
newDF <- with(DF, data.frame(YEAR = unique(YEAR)))
# prediction of smoothed estimates at each unique year value
# with standard error    
B <- predict(A,  newDF, type="response", se.fit=TRUE)


# finite difference approach to derivatives following
# example from ?predict.gam

eps <- 1e-7
X0 <- predict(A, newDF, type = 'lpmatrix')


newDFeps_p <- newDF + eps

X1 <- predict(A, newDFeps_p, type = 'lpmatrix')

# finite difference approximation of first derivative
# the design matrix
Xp <- (X0 - X1) / eps

# first derivative
fd_d1 <- Xp %*% coef(A)

# second derivative
newDFeps_m <- newDF - eps

X_1 <- predict(A, newDFeps_m, type = 'lpmatrix')
# design matrix for second derivative
Xpp <- (X1 + X_1 - 2*X0)  / eps^2
# second derivative
fd_d2 <- Xpp %*% coef(A)

If you are using boot strapping to get the confidence intervals, you should be able to get confidence intervals on these approximations.