Estimating point values in ArcGIS Pro using another point layer

You could use IDW in python (without Spatial Analyst) with the GDAL library which is included in the standard ArcGIS Pro python install (arcgispro-py3).

Here is some example code that uses gdal.Grid IDW interpolation to generate an interpolated raster from the green points, sample the interpolated values under the red points and then writes out a CSV file that you could then join back to the red points. For more information on GDAL gridding interpolation options see the gdal_grid documentation and GDAL Grid Tutorial.

from math import floor, ceil
from osgeo import (gdal, ogr)
import pandas as pd

gdal.UseExceptions()

def sample(raster, points, id_field="ID", sample_field="VALUE", layer=None):
    ras_ds=gdal.Open(raster)
    gt_forward=ras_ds.GetGeoTransform()
    gt_reverse=gdal.InvGeoTransform(gt_forward)
    rb=ras_ds.GetRasterBand(1)

    pt_ds=ogr.Open(points)

    if layer is None:
        lyr = pt_ds.GetLayer()
    else:
        lyr = pt_ds.GetLayerByName(layer)

    data = {
        id_field: [],
        sample_field: []
    }

    for feat in lyr:
        geom = feat.GetGeometryRef()
        mx,my=geom.GetX(), geom.GetY()  #coord in map units

        #Convert from map to pixel coordinates.
        px, py = gdal.ApplyGeoTransform(gt_reverse, mx, my)
        px = floor(px) #x pixel
        py = floor(py) #y pixel

        intval = rb.ReadAsArray(px, py, 1, 1)
        data[id_field].append(feat.GetField(id_field))
        data[sample_field].append(intval[0][0])  #intval is a numpy array, length=1,1 as we only asked for 1 pixel value

    return pd.DataFrame(data)


def interpolate(input_points, input_field, output_raster, pixel_size, input_layer=None):

    pt_ds = ogr.Open(input_points)
    if input_layer is None:
        lyr = pt_ds.GetLayer()
        layers = None
    else:
        lyr = pt_ds.GetLayerByName(input_layer)
        layers = [input_layer]

    xmin, xmax, ymin, ymax = lyr.GetExtent()
    ncols = ceil((xmax - xmin) / pixel_size)
    nrows = ceil((ymax - ymin) / pixel_size)

    gdal.Grid(
        output_raster,
        input_points,
        layers=layers,
        zfield=input_field,
        algorithm="invdist",
        format='GTiff',
        width=ncols,
        height=nrows
    )


if __name__ == "__main__":
    # Example using shapefiles
    input_points = "green.shp"
    input_field = "VALUE"

    sample_points = "red.shp"
    sample_field = "VALUE"
    id_field = "ID"
    output_raster = "idw.tif"

    output_csv = "values.csv"

    pixel_size = 50

    interpolate(input_points, input_field, output_raster, pixel_size)

    samples = sample(output_raster, sample_points, id_field, sample_field)
    samples.to_csv(output_csv, index=False)

    # Example using FileGDB feature classes
    input_points = "points.gdb"
    input_layer = "green"
    input_field = "VALUE"

    sample_points = "points.gdb"
    sample_layer = "red"
    sample_field = "VALUE"
    id_field = "ID"

    output_raster = "idw.tif"
    output_csv = "values.csv"

    pixel_size = 50

    interpolate(input_points, input_field, output_raster, pixel_size, input_layer=input_layer)

    samples = sample(output_raster, sample_points, id_field, sample_field, sample_layer)
    samples.to_csv(output_csv, index=False)

And here is the output (using QGIS as I'm not at work today and don't have access to ArcGIS Pro):

values.csv: