It is bit tricky, especially at raster edges so I came up with this:
select avg(value), cornerid
from corners_final cf
join lateral (select st_union(st_clip(rast, st_buffer(cf.geom, 10))) raster from dem_slope where rast && st_buffer(cf.geom, 10)) slope on true
join lateral (select unnest(st_neighborhood(raster, cf.geom, 1,1)) value) vals on true
group by cornerid
It is necessary to retrieve all rasters that are within a 3x3 window around your points, to be able to compute the neighborhood later on. This is done by firs join in the query. Please adjust the buffer amount according to your requirements, it should be a bit bigger than 3x3 pixel diameter. So after first join, each point has its own raster, that contains only pixels around them.
Then, in second join, st_neighborhood function is used to retrieve values from the desired neighborhood. Because st_neighborhood returns 2d array, unnest is necessary to retrieve individual pixel values. After that, it is only a matter of siple everage of all values for a specified point (cornerid). Group by is added only for convenience, this way you can specify multiple corner id's, or ewen avoid that condition and still be able to retrieve average values per corner id.
I have tested this query on different datasets (municipality points and DEM raster), and it performs pretty good... about 200ms/100 points.
It is also worth notice, that if you have small rasters, then it may be faster to simply union them without clipping...
Hope that helps.