Using python to create cuboid with a hole in the centre
Using from_pydata

Make  an inner and outer square of vertices for top and bottom.  Fill top  with  4 trapezium faces.   All other face indexes can be calculated from the top faces by adding appropriate index offset to them.



import bpy
from mathutils import Vector
context = bpy.context

scene = context.scene
mesh = bpy.data.meshes.new("mesh")
# 2d square
def xysquare2d(z=0, scale=1):
    r = [-scale, scale]
    return [Vector((x, y, z)) for x in r for y in r]

verts = []
# top verts
verts.extend(xysquare2d(z=2, scale=2))
verts.extend(xysquare2d(z=2))
# bottom verts       
verts.extend(xysquare2d(z=-2))
verts.extend(xysquare2d(z=-2, scale=2))
top_faces = [(0, 1, 5, 4), 
         (1, 3, 7, 5),
         (3, 2, 6, 7),
         (2, 0, 4, 6)]
# make all the other faces from top faces
bottom_faces = [(i + 8, j + 8, k + 8, l + 8) for i, j, k, l in top_faces]
outer_faces =  [(i, j, k + 8, l + 8) for i, j, k, l in top_faces]
inner_faces =  [(i + 8, j + 8, k, l) for i, j, k, l in top_faces]
# add them 
faces = top_faces  + bottom_faces + outer_faces + inner_faces

mesh.from_pydata(verts, [], faces)
obj = bpy.data.objects.new("obj", mesh)
scene.objects.link(obj)


Can create pydata from a mesh.
Make the model, then run this script on the mesh.

print("verts = [")
for v in mesh.vertices:
    print("\t", [axis for axis in v.co], ",")
print("\t]")
print("faces = [")
for f in mesh.polygons:
    print("\t", [i for i in f.vertices], ",")    
print("\t]")


Which when run on mesh created above results in 

verts = [
     [-2.0, -2.0, 2.0] ,
     [-2.0, 2.0, 2.0] ,
     [2.0, -2.0, 2.0] ,
     ....
     [2.0, 2.0, -2.0] ,
    ]
faces = [
     [0, 1, 5, 4] ,
     [1, 3, 7, 5] ,
     ....
     [10, 8, 4, 6] ,
    ]


Using boolean modifier

If you don't mind ngons, can use boolean difference, append this script to your original

# copy the mesh
grat2obj = grat1obj.copy()
grat2obj.data = grat1mesh.copy()
scene.objects.link(grat2obj)
# move it up 1, 1
grat2obj.location += Vector((-0.5, 1.25, 1.25))
# scale it 
grat2obj.scale = (2, 0.5, 0.5)

# add a boolean modifier
bm = grat1obj.modifiers.new("Chop", 'BOOLEAN')
bm.object = grat2obj
bm.operation = 'DIFFERENCE'
# apply the modifier
scene.objects.active = grat1obj
bpy.ops.object.modifier_apply(modifier=bm.name)
# remove chopper
scene.objects.unlink(grat2obj)


Using bmesh

Add a plane, inset it, remove middle face, solidify.

import bpy
import bmesh
mesh = bpy.context.object.data
bm = bmesh.new()
bmesh.ops.create_grid(bm,
                      x_segments=2,
                      y_segments=2,
                      size=1)

bmesh.ops.inset_individual(bm, faces=bm.faces, thickness=0.3)
# take out middle
[bm.faces.remove(f) for f in bm.faces if all(not e.is_boundary for e in f.edges)]                  
# extrude
bmesh.ops.solidify(bm, geom=bm.faces, thickness=1) 
bm.to_mesh(mesh)
mesh.update()