I drew a figure on GeoGebra to explore the area of the smaller quadrilateral formed by joining the vertices of the larger quadrilateral and the trisect points on the edges of the larger quadrilateral.
Here is what I think to be correct: $$Area_{\Delta E_2H_2G_2F_2}=\frac{2}{5}Area_{\Delta A_1B_1C_1D_1}$$
I haven't figured out a way to prove this question yet. If this is correct, could any of you provide a possible solution to the problem?
Let us introduce a simplified system of letters on the following figure:
Fig. 1 : See explanations below for the coordinates of $A,B,C,D$. The figure on the right illustrates a somewhat non-intuitive fact : taking a point $D$ symmetrical to its initial position (figure on the left) vs. line bissector $y=x$ gives rise to a quadrilateral $ABCD$ keeping, evidently, the same area, but to an interior quadrilateral $EFGH$ which isn't an isometric image of the initial $EFGH$, therefore with a different area in general.
Let $$r:=\dfrac{Area_{EHGF}}{Area_{ABCD}}$$
Alas, the result
$$r=\frac{2}{5}$$
is false in general (see below).
This paradox has been well underlined for a very similar case in the following reference : "Fooled by rounding" and the connected one where the "original sin" is to consider only quadrilaterals that are trapezoidal or close to be trapezoids.
In order to have a global view, let us fix 3 of the four points :
$$A=(0,1), \ \ B=(0,0), \ \ C=(1,0)$$
and let $D=(x,y)$ vary in the domain defined by $x+y \geq 1$ (in order that quadrilateral $ABCD$ is convex). The important thing is that we can do this WLOG because being given any quadrilateral, there exists an affine map that "sending" this quadrilateral to a particular quadrilateral where $ABC$ is a right isosceles triangle, with preservation of the ratio of areas (because the ratio of areas is an affine invariant (see slide 5 of this good classification).
Here is (Fig. 2) a graphical representation of function $(x,y) \mapsto r(x,y)$ as a surface with its level lines.
Fig. 2 : Representation of $r=r(x,y)$ as a surface with some of its level lines.[The tiny triangular region one finds on the left, for which $x+y<1$ isn't significant : $D$ cannot be there]. There is a "notorious" point $S=(1,1,0.4)$ featuring the case where $ABCD$ is a square (in fact all rectangles are mapped to this point as well). It is not the only point with altitude $1$ ; in fact, there is line crest of points sharing with $S$ the altitude 0.4. Another special position is $D=(0.5,0.5)$ where $r=19/48=0.3958...$ (case where $A,D,C$ aligned). The range of values of $r$ looks bounded from below with a least upper bound around $0.38$ (on an experimental basis). See Fig. 3 for an explicit calculation.
Fig. 3 : An explicit calculation of $r$ in a "coalescence" case : $D$ coincides with $A$ ; therefore $E$ and $F$ coincide with them. In this case, it is not difficult to show that $r=\frac{8}{21}\approx 0.381$ (using the fact that the abscissa of $G$ is $4/7$).
One way to verify if the proposed relationship is true.
Whether with pure metric geometry or with Cartesian coordinates of the vertices, sides and diagonals of the two quadrilaterals involved, after a certain calculation (which can be tedious) can be assumed as problem data.
Now, if $a, b, c, d$ are the sides of the greater quadrilateral, $l$ the diagonal that starts from the intersection of sides $a$ and $b$, and the other diagonal $l_1$ starting from the intersection of sides $a$ and $d$, you have the following formula for the area of the quadrilateral $$A=\dfrac12\sqrt{(2(ll_1)^2-(a^2-b^2+c^2-d^2)^2}$$ and naturally for the smaller quadrilateral
$$A'=\dfrac12\sqrt{(2(l'l_1')^2-(a_1^2-b_1^2+c_1^2-d_1^2)^2}$$
Thus we are done if we verify that $$(2(l'l_1')^2-(a_1^2-b_1^2+c_1^2-d_1^2)^2=\dfrac{4}{25}\left((2(ll_1)^2-(a^2-b^2+c^2-d^2)^2\right)$$
All the magnitudes in this formula are easily calculable but somewhat tediously.
$$XXXXXXXXXX$$
COMMENT.-$(1)$ The coordinates of $A_1,B_1,C_1, D_1$ determine the involved trisected points $A'_1,B'_1,C'_1, D'_1$.
$(2)$ Lines $A_1A'_1,B_1B'_1,C_1C'_1,D_1D'_1$ being calculated, by intersection of right lines we can determine the new vertices $E_2,F_2,G_2,H_2$.
$(3)$ The requested area, among other ways, can be calculated the following way: choosing any vertex, say $E_2=(x_E,y_E)$ and the other vertices anti-clockwise we form the following table $$x_E\hspace{10mm}y_E\\x_F\hspace{10mm}y_F\\x_G\hspace{10mm}y_G\\x_H\hspace{10mm}y_H\\x_E\hspace{10mm}y_E$$ where there are four products in descended arrows $(x_Ey_F+x_Fy_G+x_Gy_H+x_Hy_E)$ and four products in ascendent arrows $(x_Ey_H+x_Hy_G+x_Gy_F+x_Fy_E)$
$(4)$ The area $A$ is given by $$A=\frac12\big((x_Ey_F+x_Fy_G+x_Gy_H+x_Hy_E)-(x_Ey_H+x_Hy_G+x_Gy_F+x_Fy_E)\big)$$
The relation given by the OP seems to be true. In the attached figure, we apply the methode above and we find out the area of the smaller quadrilateral equal to $21.5$ which is not equal to the $\dfrac25$ of the area of the larger quadrilateral which gives $21.16$. However, the small difference should be due to the approximate decimals of the coordinates in play.
Not a complete answer, but an approach
A shearing transformation (e.g., $(x, y) \mapsto (x + ay, y)$) will preserve areas and proportions along lines. You can shear in $x$ to make $A_1D_1$ vertical, and then shear in $y$ to make $A_1B_1$ horizontal.
Scaling in either $x$ or $y$ will preserve area-proportions, so you can do that to make the distances from $A_1$ to both $B_1$ and $D_1$ be $3$.
So the problem reduces to this:
Given a quadrilateral with vertices at $(0,0), (3, 0), (3x, 3y), (0, 3)$ (with $x, y > 0$), show that the quad $Q$ bounded by the segments from
has $2/5$ the area of the quad with those original four vertices, which is, uh... \begin{align} A &= 9 + 3(3y-3)/2 + 3 (3x-3)/2 \\ &= 9 + 9(y-1)/2 + 9(x-1)/2 \\ &= 9(1 + (y-1)/2 + (x-1)/2) \\ &= 9(1 + y/2-1/2 + x/2-1/2) \\ &= 9((x+y)/2)) \end{align} Double-checking: if $x = y = 1$, so that it's a rectangle, that gives exactly $9$. Good.
So you want to show that the area of the small quad is $$ A' = \frac{2}{5}\frac{9}{2} (x+y) = \frac{9}{5} (x+y). $$
And you could either compute the exact intersections of the segments to find the area of the smaller quad, or you could apply a pair of shearing transformations to it to make it axis-aligned and compute its area as I just did for the larger quad.
Some example caculation for retangular case. Let $ABCD$ be a rectangular with $AB=6,AD=3$.
Let $TQ=x$,$XM=y$ and $\angle ZTX=\alpha$, then $XM=ZP=x$; $TQ=YN=y$; $DX=BZ=3x$;$AT=CY=3y$ and $XYZT$ is parallelogram rhombus with $XY=6x$ and $XT=6y$
We have $37a^2=(6a)^2+a^2=AB^2+AQ^2=BQ^2=100x^2$, then $x=\frac{\sqrt{37}}{10}a$ $13a^2=(2a)^2+(3a)^2=AD^2+DM^2=AM^2=(10y)^2=100y^2$, then $y=\frac{\sqrt{13}}{10}a$
and $cos(\alpha)=\frac{TQ^2+TA^2-AQ^2}{2TQ.TA}=\frac{x^2+9y^2-a^2}{6xy}=\frac{9}{\sqrt{481}}$, so $sin(\alpha)=\frac{20}{\sqrt{481}}$
Then $S_{XYZT}=sin(\alpha)XY.XT=\frac{20}{\sqrt{481}}.36\frac{\sqrt{481}a^2}{100}=7.2a^2$, and $S_{ABCD}=18a^2$ so $\frac{S_{XYZT}}{S_{ABCD}}=\frac{2}{5}$
