I'm not sure if this is the best option, but in terms of graphics it would be interesting to plot and compare both continuous and discrete PDF's and CDF's, as well as contour plots.
Let me show what I usually do:
randomWalk[x_] := RandomVariate[NormalDistribution[0, 1], x]
index = Table[randomWalk[10000], {2}];
asset[i_] := index[[i]];
pair1 = Transpose[{asset[1], asset[2]}];
Now the discrete bivariate histogram:
G1 = Histogram3D[pair1, {0.25}, "PDF", ColorFunction -> "Rainbow"]

Now the continuous bivariate plot:
G2 = Plot3D[Evaluate@PDF[BinormalDistribution[0/1],{x, y}],{x,-3,3},{y,-3,3},ColorFunction->"Rainbow"]

The countour plot:
G3 = ContourPlot[PDF[BinormalDistribution[0/1],{x,y}],{x,-3,3},{y,-3,3},ColorFunction->"Rainbow"]

Plotting the CDF's. The discrete one:
G4 = Histogram3D[pair1, {0.25}, "CDF", ColorFunction -> "Rainbow"]

The continuous one:
G5 = Plot3D[CDF[BinormalDistribution[1/2],{x, y}],{x,-4,4},{y,-4, 4},ColorFunction->"Rainbow"]

And the CDF's contour plot:
G6 = ContourPlot[CDF[BinormalDistribution[1/2],{x,y}],{x,-4,4},{y,-4,4},ColorFunction->"Rainbow"]

P.S.: you can even mix the plots, showing them at the same time for comparison. You can try, for instance
Show[G1, G2]

Show[G4, G5]

UPDATE
Using your functions to plot the PDF, CDF and contour plots. Consider:
pdf[x_, y_] := x^2 y E^(-x (y + 1)) UnitStep[x] UnitStep[y];
cdf[x_, y_] := (1 - E^-x + ((y x + 1) E^(-x (y + 1)))/(y + 1) + (y E^(-x (y + 1)) - 2 y - 1)/(y + 1)^2) UnitStep[x] UnitStep[y];
Now the plots:
Plot3D[Evaluate@pdf[x, y], {x, 0, 4}, {y, 0, 6}, ColorFunction -> "Rainbow", PlotRange -> All, Mesh -> Full, Exclusions -> None]

The PDF's contour plot:
ContourPlot[pdf[x, y], {x, 0, 4}, {y, 0, 6}, ColorFunction -> "Rainbow"]

The CDF:
Plot3D[Evaluate@cdf[x, y], {x, 0, 7}, {y, 0, 15},ColorFunction -> "Rainbow", PlotRange -> All, Mesh -> Full, Exclusions -> None]

And, finally, the CDF's contour plot:
ContourPlot[cdf[x, y], {x, 0, 7}, {y, 0, 15}, ColorFunction -> "Rainbow"]

I hope this is useful.