After thinking about it I must say it is not as simple as I thought it would be. The JW transformation on the transverse Ising model contains quite a few subtleties.
So to proceed,
1) Take your ground state for ANY $h$ expressed in the spinless fermion language. I stress ANY because this condition is true always - it's not just for $h<1$. Now this is the vacuum, specified by $| 0 \rangle$ s.t. $a_k |0\rangle = 0$. This is a non-trivial condition written in the spin-language, i.e. we take the operators $a_k$, and do the following:
2) Apply the reverse Bogoliubov transform on $a_k$: $\{a_k\}\to \{b_k\}$.
3) Apply an inverse Fourier transform: $\{b_k\} \to \{b_i\}$
4) Apply the inverse Jordan Wigner transformation: $b_i = f(\sigma^x,\sigma^y,\sigma^z)$ .
All these transformations are invertible (see this pdf for JW and inverse JW transformation), which is why you can do that. So composing all the maps one can express $a_k = g(\sigma^x, \sigma^y, \sigma^z)$, where $g$ is the highly non-trivial function.
Then one has to find the kernel of $g(\sigma^x, \sigma^y, \sigma^z)$, i.e. $|\psi\rangle $ s.t. $g(\sigma^x, \sigma^y, \sigma^z) |\psi\rangle = 0$. $|\psi \rangle $ is the ground state written in the spin-basis.
You can write a program to do it for you symbolically, but for all your effort what you'll end up with is a highly non-local ground state in the spin-basis because of all the Jordan-Wigner strings.
Remark:
There are many subtleties associated with this transformation. What is very often not mentioned when one derives the spectrum $\epsilon(k,h)$ is that the JW transformation MUST be performed separately on states with different parity in the Hilbert space.
This is because the imposition of periodic boundary conditions in spin space implies the imposition of periodic boundary conditions for ODD number of fermions but anti-periodic boundary conditions for EVEN number of fermions. This affects the Fourier transform. In the calculation of any macroscopic quantity in the thermodynamic limit, there is no difference, and many books/resources just discard talking about the two cases. But this distinction must be made if one wants to be careful about it.
One question that arose to me when I was thinking about this problem was: hm, in one limit, $h\to\infty$, the ground states is unique, while in the other limit $h \to 0$ the ground states is 2-fold degenerate. Can I see that in the fermion language easily?
There are two resolutions I can think to the problem.
1) It could be that the ground state $|0\rangle_k$ for each $k$ is not unique. That is, instead of the irreducible 2-dimension representation of the fermionic CAR that we usually assume $a_k$ acts in, $a_k$ could be operators in a $2 \times d$ (reducible) dimensional representation, with $d$ 'ground states'.
2) The even and odd sectors of the full Hilbert space give rise to two conditions on the ground state: $a_k$ in the even sector gives a condition $g(\sigma^x, \sigma^y, \sigma^z) |\psi \rangle = 0$ while $a_k$ in the odd sector gives another condition $g'(\sigma^x, \sigma^y, \sigma^z)|\psi'\rangle$ = 0.
It could perhaps be the case that when $h\to \infty$, $|\psi\rangle = |\psi'\rangle$, while when $h \to 0$, $|\psi\rangle \neq |\psi'\rangle$.
It would seem more likely to me that 2) is the correct analysis, though it's going to be one tough assertion to prove.
Basically, the answer is yes: $H$ is TRI because it is real. Reality condition really means that the Hamiltonian obeys a certain anti-unitary symmetry. In this case, the time-reversal operation is simply $T=K$ where $K$ is the complex conjugation. It is not the usual one($T=K\prod_i i\sigma^y_i$), and in particular $T^2=1$, so there is no Kramers' theorem and the spectrum is not doubly degenerate. The fact that level statistics follows GOE of course is a consequence of the reality condition. In fact, I think if there was a $T^2=-1$ time-reversal symmetry, the statistics would follow a different ensemble (the symplectic one, I believe).
You asked what if one changes the transverse field $g\sum_i \sigma^x_i$ to the $y$ direction. In that case, the two Hamiltonians are unitarily related (i.e. a $\pi/2$ spin rotation $U_z$ around $z$ would bring it back). Let me call the Hamiltonian with $x$ transverse field $H_x(g)$ where $g$ is the transverse field, and with $y$ transverse field $H_y(g)$. Define $U_z=e^{i\pi \sum_i\sigma^z_i/4}$, then it is easy to check that $U_z H_y(g) U_z^{-1}= H_x(g)$. Since we know $H_x^*(g)=H_x(g)$, we can easily find $H_y^*=U_z^2 H_y U_z^{-2}$. Therefore one just has to redefine the time-reversal symmetry to be $T=K U_z^2$. If you really want to break the reality condition, in a way that can not be fixed by additional unitary transformations, then one needs to turn on transverse fields along all three dimensions.
Last comment on your "Arguments for yes": the first argument you gave, namely one also flips the external parameters, does not work. In this way, there would be no time-reversal symmetry breaking, except the CP violation in the fundamental processes! When we talk about the symmetry of a Hamiltonian, we should just treat the system on its own, not with all the external devices that generate the various terms -- unless you want to consider the dynamics of these devices, but then it is a different Hamiltonian.
Best Answer
This is a very good question. The same operator algebra does not imply that $H(J,h)$ and $H(h,J)$ have the same spectrum. As has been mentioned in Dominic's answer, even the ground state degeneracy is different under the interchange of $J$ and $h$ ($J\gg h$: symmetry-broken two-fold degeneracy, and $J\ll h$ unique ground state), therefore it is impossible to establish a one-to-one mapping between the eigenstates of $H(J,h)$ and $H(h,J)$. One must keep in mind that the duality transformation only preserves local dynamics but the global (topological) properties will be lost. This statement becomes sharper in higher dimensions. Like in 2D, the $\mathbb{Z}_2$ lattice gauge theory is dual to the quantum Ising model, however the topological order of the gauge theory (most prominently the topology-dependent ground state degeneracy) is completely lost in the dual Ising model, even though there is a beautiful correspondence between their local excitations (e.g. charge and vison).
Nevertheless, duality is still very useful if we only focus on the local excitations. Many important problems like low-energy dynamics, phase transitions and criticality are only related to local excitations, then the duality transformation can help us a lot in understanding these things.
To appreciate the duality in the 1D transverse field Ising model, it is better to just look locally and try to figure out the local correspondence of the bulk excitations, without caring too much about the global properties such as boundary conditions, infinite strings, ground state degeneracy etc. The idea of duality is actually simple: one can describe a 1D Ising chain either by the spin variables living on each site, or by the kink variables living on each link. A kink on the spin chain is a link across which the Ising spins are opposite. Then every link $l$ can only have two possible states: $$\tau_l^z=\left\{\begin{array}{cc}+1 & \text{unkinked,}\\ -1 & \text{kinked,}\end{array}\right.$$
If we have specified all the kink configuration $\tau_l^z$ on each link $l$, we can actually determine the spin configuration $\sigma_i^z$ on each site $i$, with only one additional piece of knowledge about the left-most spin $\sigma_0^z$. The trick is just to accumulate the kink configurations from left all the way to the right, $$\sigma_i^z=\sigma_0^z\prod_{0<l<i}\tau_l^z.$$ So the spin configuration is uniquely determined by the kink configuration (up to the left-most spin). If we count the number of states in the Hilbert space, the Hilbert space dimension will be $2^{N_\text{site}}$ in the spin language and $2^{N_\text{link}}$ in the kink language, where $N_\text{site}$ and $N_\text{link}$ are the number of sites and links respectively, which are equal (apart from the left-most site) on the 1D lattice, so the Hilbert space dimension is actually the same in both languages. In this sense, we can say that the correspondence between the spin and the kink descriptions is almost one-to-one (especially in the thermodynamic limit), even though there might be some complication arising from the left-most boundary (which however is to be ignored in the duality transform, as duality only cares local properties).
Now we can rephrase the original transverse field Ising model $$H=-J\sum_{i}\sigma_i^z\sigma_{i+1}^z-h\sum_{i}\sigma_i^x,$$ in the kink language. It is not hard to see that the coupling between the Ising spins is just the chemical potential for the kink $$-J\sum_{i}\sigma_i^z\sigma_{i+1}^z=-J\sum_l\tau_l^z,$$ which basically follows from the physical meaning of the kink variable $\tau_l^z$. Admittedly this equality may run into a little trouble on the boundary, where some sites or links might be missing, but they do not affect the local properties deep in the bulk, so we just ignore. The translation of transverse field term is more involved. In the spin language $\sigma_i^x$ operator just flips the spin on site $i$, which would correspond to simultaneous creation or annihilation two kinks on the links adjacent to that site, or moving an existed kink across the site.
In either case, flipping a spin would correspond to simultaneously changing the kink variables $\tau^z$ on adjacent links, which can be carried out by $\tau^x_l\tau^x_{l+1}$, s.t. $$-h\sum_i\sigma_i^x=-h\sum_{l}\tau^x_l\tau^x_{l+1}.$$ Obviously the relation between $\tau^x$ and $\tau^z$ is exactly the same as our familiar $2\times 2$ Pauli matrices $\sigma^x$ and $\sigma^z$, s.t. $\tau^x|\tau^z=+1\rangle=|\tau^z=-1\rangle$ and $\tau^x|\tau^z=-1\rangle=|\tau^z=+1\rangle$, and the algebraic relations like $\{\tau^x,\tau^z\}=0$ are just consequences that follow. As the OP have pointed out, the algebraic relations can not guarantee the representation to be fundamental, it is actually the above physical picture that guarantees the representation of the kink operators.
Putting together the above results, we arrive at the Hamiltonian in terms of the kink operators $\tau_l^x$ and $\tau_l^z$ $$H=-h\sum_{l}\tau^x_l\tau^x_{l+1}-J\sum_l\tau_l^z.$$ One might have been wondering for a while that why I kept using the symbol $\tau$ instead of $\mu$ in the original post. Now it is clear that $\tau$ is still one step away from $\mu$ by a basis transformation on each link that redefines $\tau^x=\mu^z$ and $\tau^z=\mu^x$ (relabeling $x\leftrightarrow z$). Such a unitary transform will not change any physics, but just to bring back the standard form of the transverse field Ising model to accomplish the duality transform, $$H=-h\sum_{l}\mu^z_l\mu^z_{l+1}-J\sum_l\mu_l^x.$$ So the duality between $J$ and $h$ is now manifest, but what is the physical meaning of $\mu_l^z$ indeed? To answer this question, one should first understand that the relation between $\tau^z$ and $\tau^x$ is just like that between the coordinate and the momentum. They are related by a $\mathbb{Z}_2$ version of the Fourier transformation. If we treat the two states of $|\tau^z=\pm 1\rangle$ as two position eigenstates in a two-site system, then the $\tau^x$ eigenstates $|\tau^z=+1\rangle\pm|\tau^z=-1\rangle$ are nothing but the momentum eigenstates with the momentum = $0$ and $\pi$ respectively. In this sense, we can say $\mu_l^z\equiv\tau_l^x$ is the conjugate momentum of the kink variable $\tau_l^z$ on each link. In fact, this concept is so important that people invent a name for $\mu^z$, i.e. the vison variable, s.t. $$\mu_l^z=\left\{\begin{array}{cc}+1 & \text{vison off,}\\ -1 & \text{vison on.}\end{array}\right.$$ By saying that the vison is the conjugate momentum of the kink, we mean that if there is a vison sitting on a link then the kinked and the unkinked configurations across that link will be differed by a minus sign in the wave function. Unlike the kink which is just another way to encode the spin configuration, the vison does not have a correspondent spin configuration. In fact, the vison configuration is encoded in the relative sign between different spin configurations in the wave function. It represents the inter-relation among the spin configurations other than any particular spin configuration itself, or in other words, the quantum entanglement in the spin wave function.
Mathematically this can be seen from the fact that the vison operator $\mu_l^z$ is non-local in terms of the spin operator $$\mu_l^z=\prod_{i<l}\sigma_i^x=\tau_l^x,$$ which flips all the spins to its left to create (or annihilate) a kink. Now let us discuss more about this infinite string of $\sigma^x$ stretching all the way to the left. A first question is that can we "gauge" this string to right? This can be done by applying the operator $S\equiv\prod_i\sigma_i^x$, as $$\mu_l^z\to S\mu_l^z= \prod_{i}\sigma_i^x\prod_{i<l}\sigma_i^x = \prod_{i>l}\sigma_i^x.$$ One can see the operator $S$ simply flips all the spins in the system, meaning that it actually implement the global Ising symmetry transformation to the spins $\sigma_i^z\to-\sigma_i^z$. Because $S$ is a symmetry of the Hamiltonian (as $[S,H]=0$), the eigenstates of $H$ are spitted into the even ($S=+1$) and the odd ($S=-1$) sectors. In the even sector, we can gauge the string to the right; while in the odd sector, gauging the string to the right will induce a $\mathbb{Z}_2$ gauge transformation of the vison $\mu_l^z\to-\mu_l^z$. So combined with the gauge transformation of the vison, the vison string can be made invisible indeed. Then what is the significance of the this vison string? Recall that $\mu_l^z=\tau_l^x$ is also the creation/annihilation operator of the the kink. So applying a string of $\sigma^x$ operators will actually create two kinks at both ends of the string $$\tau_{l_1}^x\tau_{l_2}^x=\prod_{l_1<i<l_2}\sigma_i^x.$$ The kink is a local excitation of the system (in the $J>h$ phase), so it can be viewed as a particle. From this perspective, we can see emergent particles at the ends of the string, which is exactly one of the central theme of Levin-Wen model and string-net condensation.
In fact, there is a very interesting relation between the 1D transverse field Ising model and the 2D Levin-Wen model (with $\mathbb{Z}_2$ topological order), that the former can be considered as a synthetic dislocation of the latter by anyon condensation, which was described in a paper (arXiv:1208.4109) that I wrote with my friend Chao-Ming and Prof. Wen. We basically showed that the 1D transverse field Ising model can emerge in the 2D $\mathbb{Z}_2$ topological ordered system as some kind of line defect, where some particular type of the strings between the anyons in 2D will naturally degrade to the vison strings along the emergent 1D Ising chain. So in this sense, the invisible string in the Ising model is really the same invisible string in the string-net condensate (but just confined to the 1D system).