I think what fundamentally needs to be explained here is this:

The physical interpretation of the Hubble time is that it gives the time for the Universe to run backwards to the Big Bang **if the expansion rate (the Hubble "constant") were constant**. Thus, it is a measure of the age of the Universe. The Hubble "constant" actually isn't constant, so the Hubble time is really only a rough estimate of the age of the Universe.

(source, emphasis added) You can verify this mathematically: if the Hubble time $1/H$ really did track the age of the universe (ignoring the general relativistic complications of what "the age of the universe" really means), then it must be the case that $H(t) = 1/t$. Given the definition of the Hubble parameter as $\dot{a}(t)/a(t)$, you can write

$$\frac{\dot{a}(t)}{a(t)} = \frac{1}{t}$$

This differential equation for $a$ has the solution $a(t) = ct$, which indicates that the universe would be expanding linearly in this case.

In reality, of course, the universe does not expand linearly, at least not always. But the available evidence suggests that it's been expanding pretty *nearly* linearly for a long time, $\dot{a}(t) \approx \text{const.}$ for the last 10-12 billion years, which is why the Hubble time is so close to the age of the universe as estimated from other methods (well, method - WMAP data).

The problem is that your calculation has no real physical meaning. It is only meaningful to compare two quantities within the same inertial frame. But there is no global inertial frame that connects us to a co-moving galaxy at the Hubble radius: such a galaxy is at rest in its own local cosmological inertial frame (ignoring local peculiar motions). So, due to the Cosmological Principle, its own cosmic proper time is similar to ours.

To obtain meaningful information, we need to analyse the light of that galaxy as we observe it. And as you indicate in your own answer, the Hubble radius is not a cosmological horizon in terms of light paths.

So see why, let us look at the Standard ΛCDM-model in more detail. The expansion of the universe can be expressed in terms of the scale factor $a(t)$ and its derivatives, with $a=1$ today. From the Friedmann equations, it can be shown that the **Hubble parameter** $H$ has the form
$$
H(a) = \frac{\dot{a}}{a} = H_0\sqrt{\Omega_{R,0}\,a^{-4} + \Omega_{M,0}\,a^{-3} + \Omega_{K,0}\,a^{-2} + \Omega_{\Lambda,0}},
$$
with $H_0$ the present-day Hubble constant, $\Omega_{R,0}, \Omega_{M,0}, \Omega_{\Lambda,0}$ the relative present-day radiation, matter and dark energy density, and $\Omega_{K,0}=1-\Omega_{R,0}- \Omega_{M,0}- \Omega_{\Lambda,0}$. I will assume the values
$$
\begin{gather}
H_0 = 67.3\;\text{km}\,\text{s}^{-1}\text{Mpc}^{-1},\\
\Omega_{R,0} = 0, \quad\Omega_{M,0} = 0.315, \quad\Omega_{\Lambda,0} = 0.685,\quad\Omega_{K,0} = 0.
\end{gather}
$$
It's important to note that, even though the expansion of the universe is currently accelerating ($\ddot{a}>0$), the Hubble parameter $H$ is in fact always *decreasing*. The cosmic time can then be calculated from
$$
\text{d}t = \frac{\text{d}a}{\dot{a}} = \frac{\text{d}a}{aH(a)},
$$
so that
$$
t(a) = \int_0^a\frac{\text{d}a}{aH(a)}.
$$
Likewise, a photon travels on a null-geodesic
$$
0 = c^2\text{d}t^2 - a^2(t)\text{d}\ell^2,
$$
with $\text{d}\ell$ a co-moving displacement, so that the **co-moving distance** travelled by a photon is
$$
D_\text{c} = c\int_{a_\text{em}}^{a_\text{ob}}\frac{\text{d}a}{a\dot{a}} = c\int_{a_\text{em}}^{a_\text{ob}}\frac{\text{d}a}{a^2H(a)},
$$
with $a_\text{em} = a(t_\text{em})$ and $a_\text{ob} = a(t_\text{ob})$ the scale-factors at the moments of emission and observation. At any moment $t$, a co-moving distance $D_\text{c}(t)$ can be converted into a **proper distance** $D(t) = a(t) D_\text{c}(t)$.

The distance to the current Hubble radius is
$$
D_\text{H}=\frac{c}{H_0}\approx 14.5\;\text{Gly}.
$$
Now, we can define two important cosmological horizons: the first is the **particle horizon**, which marks the edge of the observable universe. It is the maximum distance that light has been able to travel to us between $t=0$ and the present day, i.e. it is as far as we can see:
$$
D_\text{ph} = c\int_0^1\frac{\text{d}a}{a^2H(a)}.
$$
Since $H(a) < a^{-2}H_0\ $ for $a<1$, we get
$$
D_\text{ph} > \frac{c}{H_0}\int_0^1\text{d}a = D_\text{H},
$$
it follows that the Hubble radius is smaller than the particle horizon, and thus part of the observable universe. In fact, $D_\text{ph}\approx 46.2\;\text{Gly}$.

The second horizon is the **event horizon**. It is the region of space from which a photon that is emitted 'today' will still be able to reach us at some point in the future. Thus
$$
D_\text{eh} = c\int_1^\infty\frac{\text{d}a}{a^2H(a)}.
$$
Since $H(a) < H_0\ $ for $a>1$, we get
$$
D_\text{eh} > \frac{c}{H_0}\int_1^\infty\frac{\text{d}a}{a^2} = D_\text{H},
$$
so the Hubble radius is also smaller than the event horizon: a galaxy on the Hubble radius can still send signals to us (or we to them). $D_\text{eh}\approx 16.7\;\text{Gly}$.

The situation is displayed in these graphs (click on 'view image' for a larger version): the first is with co-moving coordinates, the second with proper coordinates.

The black solid lines indicate our present position. The blue line is the particle horizon through time, the red line is the event horizon, the green area is the Hubble sphere. The dotted black line is a co-moving galaxy that is currently located on the Hubble radius. Its photons that we observe today have travelled on the cyan path (they were emitted at $t=4.3\;\text{Gy}$). Its photons that it emits today ($t=13.8\;\text{Gy}$) will travel on the purple path (they will reach us at $t=49\;\text{Gy}$).

So can we say anything about **time dilation**? The answer is yes. The light that we observe is also redshifted
$$
1 + z = \frac{\lambda_\text{ob}}{\lambda_\text{em}} = \frac{\nu_\text{em}}{\nu_\text{ob}}.
$$
Within a small time $\delta t_\text{em}$ the source emits a light wave with $\nu_\text{em}\delta
t_\text{em}$ oscillations. Those same oscillations are observed within time $\delta t_\text{ob}$ with frequency $\nu_\text{ob}$, in other words $\nu_\text{em}\delta
t_\text{em} = \nu_\text{ob}\delta t_\text{ob}$, so that
$$
1 + z = \frac{\delta t_\text{ob}}{\delta t_\text{em}}.
$$
In other words, cosmic redshift is directly related to time dilation. Also, within a small time $\delta t_\text{em}$, the distance that light has to travel doesn't change:
$$
\int_{t_\text{em}}^{t_\text{ob}}\frac{c\,\text{d} t}{a(t)} =
\int_{t_\text{em} + \delta t_\text{em}}^{t_\text{ob} + \delta t_\text{ob}}\frac{c\,\text{d} t}{a(t)},
$$
or
$$
\int_{t_\text{ob}}^{t_\text{ob} + \delta t_\text{ob}}\frac{c\,\text{d} t}{a(t)} =
\int_{t_\text{em}}^{t_\text{em} + \delta t_\text{em}}\frac{c\,\text{d} t}{a(t)}.
$$
in these small intervals, the integrands remain constant, so that
$$
\frac{\delta t_\text{ob}}{a(t_\text{ob})} = \frac{\delta t_\text{em}}{a(t_\text{em})}
$$
and
$$
1 + z = \frac{a(t_\text{ob})}{a(t_\text{em})}.
$$
The light from the co-moving galaxy at the current Hubble radius that we observe today was emitted when $a(t_\text{em})=0.403$, so that $z=1.48$, and the observed events from this galaxy are time-dilated by a factor $2.48$. And the light it emits today will be observed when $a(t_\text{ob})=8.07$, with redshift $z=7.07$.

Such time dilations have indeed been observed: we literally see distant supernovae explode in slow motion!

## Best Answer

The general solution works as follows:

We start with the Friedmann equation $$ \dot{a}^2 - \frac{8\pi G}{3}\rho a^2 = -kc^2, $$ with $k=0,\ 1,\ $or $-1$, and $\rho$ the total density. Since the right-hand side is constant, we can write $$ \dot{a}^2 - \frac{8\pi G}{3}\rho a^2 = \dot{a}_0^2 - \frac{8\pi G}{3}\rho_0 a_0^2, $$ where the subscript 0 denotes the present-day values. If we introduce the Hubble constant $$ H_0 = \frac{\dot{a}_0}{a_0} $$ and the present-day critical density $$ \rho_{c,0} = \frac{3H_0^2}{8\pi G}, $$ we get $$ \frac{\dot{a}^2}{a_0^2} - H_0^2\frac{\rho}{\rho_{c,0}} \frac{a^2}{a_0^2} = H_0^2 - H_0^2\frac{\rho_0}{\rho_{c,0}} $$ or $$ H^2 = \frac{\dot{a}^2}{a^2} = H_0^2\left[\frac{\rho}{\rho_{c,0}} + \frac{a_0^2}{a^2}\left(1 - \frac{\rho_0}{\rho_{c,0}}\right)\right]. $$ Now, there are three contributions to the total density: radiation, matter (normal and dark) and dark energy: $$ \rho = \rho_R + \rho_M + \rho_{\Lambda}. $$ These densities change over time as follows: the matter density decreases as the volume of the universe increases, so $\rho_M\sim a^{-3}$, as you'd expect. The radiation falls off as $\rho_R\sim a^{-4}$ (the extra factor is due to redshift). And in the Standard Model, the dark energy remains constant: $\rho_{\Lambda} = \text{const}$. In other words, $$ \begin{align} \rho_R a^4 &= \rho_{R,0}\, a_0^4,\\ \rho_M a^3 &= \rho_{M,0}\, a_0^3,\\ \rho_\Lambda &= \rho_{\Lambda,0}, \end{align} $$ and finally, with the notations $$ \Omega_{R,0} = \frac{\rho_{R,0}}{\rho_{c,0}},\quad \Omega_{M,0} = \frac{\rho_{M,0}}{\rho_{c,0}},\quad \Omega_{\Lambda,0} = \frac{\rho_{\Lambda,0}}{\rho_{c,0}},\\ \Omega_{K,0} = 1 - \Omega_{R,0} - \Omega_{M,0} - \Omega_{\Lambda,0}, $$ we find $$ H(a) = H_0\sqrt{\Omega_{R,0}\,a^{-4} + \Omega_{M,0}\,a^{-3} + \Omega_{K,0}\,a^{-2} + \Omega_{\Lambda,0}}, $$ where we used the convention $a_0=1$. Also note that $$ \dot{a} = H_0\sqrt{\Omega_{R,0}\,a^{-2} + \Omega_{M,0}\,a^{-1} + \Omega_{K,0} + \Omega_{\Lambda,0}\,a^2},\\ \ddot{a} = -\frac{1}{2}H_0^2\left(2\,\Omega_{R,0}\,a^{-3}+\Omega_{M,0}\,a^{-2} -2\,\Omega_{\Lambda,0}\,a\right). $$ The latest values of the parameters, obtained from the Planck mission, are $$ H_0 = 67.3\;\text{km}\,\text{s}^{-1}\text{Mpc}^{-1},\\ \Omega_{R,0} = 9.24\times 10^{-5},\qquad\Omega_{M,0} = 0.315,\\ \Omega_{\Lambda,0} = 0.685,\qquad\Omega_{K,0} = 0. $$ So now we have the Hubble parameter as a function of the scale radius $a$. How can we convert this into a function of time? From $$ \dot{a} = \frac{\text{d}a}{\text{d}t} $$ we get $$ \text{d}t = \frac{\text{d}a}{\dot{a}} = \frac{\text{d}a}{aH(a)} = \frac{a\,\text{d}a}{a^2H(a)}, $$ so that $$ \begin{align} t(a) &= \int_0^a \frac{a'\,\text{d}a'}{a'^2H(a')}\\ &= \frac{1}{H_0}\int_0^a \frac{a'\,\text{d}a'}{\sqrt{\Omega_{R,0} + \Omega_{M,0}\,a' + \Omega_{K,0}\,a'^2 + \Omega_{\Lambda,0}\,a'^4}}. \end{align} $$ Inverting this relation, we get $a(t)$. Unfortunately, this inversion has to be done numerically. And finally, $$ H(t) = H(a(t)). $$

P.S. The solution that you mentioned is the case where radiation and matter are negligible, and dark energy has a more general form (called quintessence): $$ \rho_R=\rho_M=0,\quad \rho_\Lambda = \rho_{\Lambda,0}\,a^{-3(1+w)}, $$ where $w=-1$ corresponds with the normal case of a cosmological constant. In this case, for a universe with no curvature, $$ H^2 = H_0^2\,a^{-3(1+w)},\qquad t(a) = \frac{1}{H_0}\int_0^a a'^{(1+3w)/2}\,\text{d}a', $$ with solution $a\sim t^{2/(3+3w)}$, for $w>-1$. Solutions with $w\leqslant-1$ have no big bang, i.e. the lower bound in the integral $t(a)$ cannot be zero.

In any case, these are not accurate descriptions of our universe, since they ignore the contributions of matter and radiation.