I ended up using an approach based on the snowfall package. It is quite simple, works really good and the point extraction function is fast as the number of cores that you can use. The approach I used was inspired by this post, and here is my reproducible example:

library(raster)
library(snowfall)

# Create date sequence
idx <- seq(as.Date("2010/1/1"), as.Date("2099/12/31"), by = "day")

# Create raster stack and assign dates
# WARNING: raster stack will be ~ 400 MB in size
r <- raster(ncol=5, nrow=5)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)))))
s <- setZ(s, idx)

# Create random spatial points
pts <- SpatialPoints(cbind(x=runif(655, -180, 180),
                           y=runif(655, -90, 90)),
                     proj4string=CRS(projection(s)))

# Extract values to a data frame - multicore approach
# First, convert raster stack to list of single raster layers
s.list <- unstack(s)
names(s.list) <- names(s)

# Now, create a R cluster using all the machine cores minus one
sfInit(parallel=TRUE, cpus=parallel:::detectCores()-1)

# Load the required packages inside the cluster
sfLibrary(raster)
sfLibrary(sp)

# Run parallelized 'extract' function and stop cluster
e.df <- sfSapply(s.list, extract, y=pts)
sfStop()

# Fix resulting data frame
DF <- data.frame(t(e.df)); row.names(DF)=NULL
DF <- cbind(getZ(s), DF)
names(DF)[1] <- "Date";
names(DF)[2:length(names(DF))] <- paste0('Point ', 1:length(pts))

# Check resulting data frame
head(DF)

I hope it can serve as a future reference for people interested in doing the same.

I ended up using an approach based on the snowfall package. It is quite simple, works really good and the point extraction function is as fast as the number of cores that you can use. The approach I used was inspired by this post, and here is my reproducible example:

library(raster)
library(snowfall)

# Create date sequence
idx <- seq(as.Date("2010/1/1"), as.Date("2099/12/31"), by = "day")

# Create raster stack and assign dates
# WARNING: raster stack will be ~ 400 MB in size
r <- raster(ncol=5, nrow=5)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)))))
s <- setZ(s, idx)

# Create random spatial points
pts <- SpatialPoints(cbind(x=runif(655, -180, 180),
                           y=runif(655, -90, 90)),
                     proj4string=CRS(projection(s)))

# Extract values to a data frame - multicore approach
# First, convert raster stack to list of single raster layers
s.list <- unstack(s)
names(s.list) <- names(s)

# Now, create a R cluster using all the machine cores minus one
sfInit(parallel=TRUE, cpus=parallel:::detectCores()-1)

# Load the required packages inside the cluster
sfLibrary(raster)
sfLibrary(sp)

# Run parallelized 'extract' function and stop cluster
e.df <- sfSapply(s.list, extract, y=pts)
sfStop()

# Fix resulting data frame
DF <- data.frame(t(e.df)); row.names(DF)=NULL
DF <- cbind(getZ(s), DF)
names(DF)[1] <- "Date";
names(DF)[2:length(names(DF))] <- paste0('Point ', 1:length(pts))

# Check resulting data frame
head(DF)

I hope it can serve as a future reference for people interested in doing the same.