You are correct. I would suspect the author of this blogpost simply didn't notice his mistake because he happened to be using only a square aspect ratio (his example code seems to use an 800×800 framebuffer). Or he is using only square aspect ratios on purpose to avoid this issue. Either way, in general, you will have to consider the aspect ratio when going from normalized device coordinates back to view space.
\$
\newcommand{\mvec}[1]{\begin{pmatrix}#1\end{pmatrix}}
\newcommand{\ivec}[1]{\left(\begin{smallmatrix}#1\end{smallmatrix}\right)}
\newcommand{\mmat}[1]{\begin{bmatrix}#1\end{bmatrix}}
\renewcommand{\vec}[1]{\mathrm{\mathbf{#1}}}
\newcommand{\mat}[1]{\mathrm{\mathbf{#1}}}
\$
Anyhow, if the goal is to generate rays in a way that matches rasterization by the hardware graphics pipeline, I would personally follow a different approach. Rather than trying to construct my rays to match the way in which I expect my projection matrix to be constructed, I would simply compute the rays from a given projection matrix \$\mat M\$. All you need is the inverse \$\mat M^{-1}\$ of that projection matrix.
Given a point on the screen in viewport (pixel) coordinates \$\ivec{x && y}^T\$, we compute the clip-space coordinates of two corresponding points, one on the near and one on the far plane:
\begin{align}
\vec p''_{near} &= \mvec{
2\cdot\frac{x}{w} - 1 \\
2\cdot\frac{y}{h} - 1 \\
-1 \\
1
} &
\vec p''_{far} &= \mvec{
2\cdot\frac{x}{w} - 1 \\
2\cdot\frac{y}{h} - 1 \\
1 \\
1
},
\end{align}
where \$w\$ and \$h\$ are the width and height of the viewport respectively. For simplicity, we'll ignore the fact that the viewport may be offset. It should be obvious how to account for such an offset in case that would be necessary.
Note that in OpenGL, pixel locations correspond to the centers of the cells of a sampling grid by default. Thus, the viewport coordinates of the pixel at row index \$j\$ and column index \$i\$ would be \$x = i + \frac{1}{2}\$ and \$y = j + \frac{1}{2}\$. gl_FragCoord will already give you the correct coordinates unless you explicitly change its coordinate system origin via layout qualifier.
We can then transform these two points back to camera-space by first multiplying with the inverse projection matrix
\begin{align}
\vec p'_{near} &= \mat M^{-1} \vec p''_{near} &
\vec p'_{far} &= \mat M^{-1} \vec p''_{far}
\end{align}
and then dehomogenizing to get euclidean coordinates:
\begin{align}
\vec p_{near} &= \frac{1}{p'_{near_w}} \cdot \mvec{
p'_{near_x} \\
p'_{near_y} \\
p'_{near_z}
} &
\vec p_{far} &= \frac{1}{p'_{far_w}} \cdot \mvec{
p'_{far_x} \\
p'_{far_y} \\
p'_{far_z}
}.
\end{align}
Knowing a point on the near and one on the far plane, we can construct a view ray representing all points which project to our pixel, e.g.,
\begin{equation}
\vec r(t) = \vec p_{near} + t \cdot (\vec p_{far} - \vec p_{near}).
\end{equation}
Note how this automatically also allows us to account for clipping on the near and far plane as everything behind the near plane will be behind the ray start (\$t < 0\$) and everything beyond the far plane will correspond to \$t > 1\$.
\$\mat M\$ can also be a combined View-Projection or Model-View-Projection matrix to directly get the ray in world- or object-space instead.
There are many ways to construct a projection matrix. And perspective projection is not the only kind of projection. Above approach will automatically and exactly match anything you can throw at it in the form of a matrix. While it might be slightly more expensive than a specialized construction for a particular kind of projection matrix, I doubt you will be able to notice an impact on performance in practice, as ray setup is unlikely to be a bottleneck for any kind of raytracing that does something close to meaningful…
