Quick keyframe frame insertion for Shape Keys with `foreach_set` from data
Could be simplified to one shape key.

In the initial state of the cube, it is dimensionless (squashed down
to a flat object). In the animation, it should be clear that the Cube
starts growing from 'a flat surface'.
Also, how can I describe the scale in Z direction, when the values are
accumulative? This gets complex already.

IMO you have totally overcomplicated it, by adding a shapekey for each data point.  For this case even Using foreach methods isn't going to help particularly, ie for 4000 points of data creating as many shapekeys and keyframing 2 points per.
Am also of the opinion this could be reduced to simply keyframing object scale, a single shapekey, or solidify thickness to name but 3.

Using solidify thickness on faces defined by vert group &quot;Group&quot;
Shapekey emulating how scale works.


Add a default cube
Transform such that it is unit size, origin at base
Copy the coords
Scale it flat, make that the basis shape,
Make the other shape key unit.

Ok now have a shapekey that if zero has zero height Basis, and if 1 unit height. Changing the range allows scaling outside the (0, 1) range. As illustrated, but appears to be hardcoded to only 10
There is really no need for the shapekeys, since the result is same as scaling in z axis, So could instead create no shapekeys and simply keyframe scale.
Test script. Creates the unit cube illustrated above
import bpy
from mathutils import Matrix
bpy.ops.mesh.primitive_cube_add()
ob = bpy.context.object
me = ob.data
me.transform(
        Matrix.Translation((0, 0, .5)) @        
        Matrix.Diagonal((.5, .5, .5, 1))
        )

coords = [None] * len(me.vertices) * 3
me.vertices.foreach_get(&quot;co&quot;, coords)
me.transform(
        Matrix.Diagonal((1, 1, 0, 1))
        )
basis = ob.shape_key_add(name=&quot;Basis&quot;)   
unit = ob.shape_key_add(name=&quot;Unit&quot;)
unit.data.foreach_set(&quot;co&quot;, coords)

Some animation basics
When you add a keyframe to a shapekey, it creates an Action if there is not already one, an fcurve for the data path being animated (if there already isn't one) and inserts a 2D point for (frame, value)
It is in populating the fcurves keyframe points where foreach_set is then used.
&gt;&gt;&gt; C.object.data.shape_keys.animation_data.action
bpy.data.actions['Key.018Action']

&gt;&gt;&gt; action = C.object.data.shape_keys.animation_data.action
&gt;&gt;&gt; for fc in action.fcurves:
...     fc.data_path
...     
'key_blocks[&quot;Unit&quot;].value'


Besides, to my understanding,
foreach_set method creates and assigns all the keyfames from a data
set at once.
This is the efficient way when you have thousands of data points. But
I just could not figure out a way to use it in shape keys

Yeah?, thanks for the tip.
Knowing where the action belongs and what the data path looks like, is the key to using the methods outlined in my answers to  links posted
Fast keyframe insertion
Editing fcurve.keyframe points in FAST mode?
Example
Using data from question

Created mesh with shapekey as above,
Make the shapekey maximum (as big as it gets from data)
Have offset data such that 2000 is frame 1, and taken one from value for example sake, to get zero at frame 1.
Created an action for the &quot;Unit&quot; shapekey and populated it.

Test script.
import numpy as np
import bpy
from mathutils import Matrix

data = (
    (2000, 1),
    (2001, 2.5),    
    (2003, 3.4),    
   
    (2020, 11.9),    
    )

bpy.ops.mesh.primitive_cube_add()
ob = bpy.context.object
me = ob.data
    
d = np.array(data)
d -= (1999, 1)
frames, values = d.T

# max range of 10 ..
#unit.slider_max = values.max()


me.transform(
        Matrix.Translation((0, 0, .5)) @        
        Matrix.Diagonal((.5, .5, .5, 1))
        )
me.transform(
        Matrix.Diagonal((1, 1, values.max(), 1))
        )
        
coords = [None] * len(me.vertices) * 3
me.vertices.foreach_get(&quot;co&quot;, coords)
me.transform(
        Matrix.Diagonal((1, 1, 0, 1))
        )
basis = ob.shape_key_add(name=&quot;Basis&quot;)   
unit = ob.shape_key_add(name=&quot;Unit&quot;)
unit.data.foreach_set(&quot;co&quot;, coords)
# animation
me.shape_keys.animation_data_create()
action = bpy.data.actions.new(&quot;Grow&quot;)
action.id_root = 'KEY'
fc = action.fcurves.new(f'key_blocks[&quot;{unit.name}&quot;].value')

values /= values.max()
fc.keyframe_points.add(count=len(frames))
fc.keyframe_points.foreach_set(&quot;co&quot;, d.ravel())

me.shape_keys.animation_data.action = action