Is there a wave function for anyons? People talk about anyons a lot. 
But i have never seen an anyon wave function. 
I suspect that there is no such thing as a wave function for anyons. I mean, anyons are not generalizations of bosons or fermions. For bosons and fermions, one can have many-body hilbert spaces, but for anyons, there is no such thing. 
 A: I don't think that is the case. A useful reference is: http://journals.aps.org/rmp/abstract/10.1103/RevModPhys.80.1083.
One way to approach a theory of anyons is to start by writing down a list of particle types along with their fusion rules. Once doing this one may obtain consistency equations from solving the hexagon and pentagon equations arising from modular tensor categories. If successful in solving these equations, you have a viable anyon theory. If multiple solutions exist you have multiple theories. 
Now once doing this, we can label each state in the hilbert with a fusion tree. Hence our hilbert space is very well defined - albeit abstract. We can provide interactions on this hilbert space by creating projectors which favor nearby (in real space) anyons to fuse together in various channels. Whence doing this, in principle we could diagonalize a finite sized system and extract the wavefunctions. 
A simple example is to consider Ising anyons. These appear in the Kitaev honeycomb model and manifest them selves as Majorana zero modes in the 1 dimensional p-wave wires. (see for example http://www.sciencedirect.com/science/article/pii/S0003491605002381, and). In the case of 1D p-wave wires we can certainly write down the wavefunctions of the Majorana zero modes as they happen to be solutions of the BdG equations. 
A: For those with a similar query, it's worth reading the original paper in which anyons were proposed, written by Leinaas and Myrheim in 1977 (https://www.ifi.unicamp.br/~cabrera/teaching/referencia.pdf). It provides a perfectly consistent 1st quantised theory of anyons. I have attempted to summarise the main claims of the paper below, but I'm not a mathematician, so do let me know if I'm getting something wrong!
The Wave Function is over Configuration Space
The authors of the paper point out that a description in which symmetry constraints arise naturally can be obtained by letting the wave function be defined over the system's classical configuration space. For example, in the configuration space of two identical particles, the point $(\mathbf{x}_1, \mathbf{x}_2)$ is identified with $(\mathbf{x}_2, \mathbf{x}_1)$ - i.e. the point where the two particles are swapped. These two points are completely indistinguishable from each other.
The reason for non-trivial exchange statistics is that the points at which the locations of particles coincide, such as $(\mathbf{x}, \mathbf{x})$, are removed from the configuration space. That is, the space has punctures in it, making it not simply connected. Note that if you do not remove these points, you will obtain bosonic exchange statistics.
The wave function is a multi-valued function from this punctured configuration space into the complex plane (or into a higher-dimensional vector space over the complex field, if you are dealing with spin and other degeneracies). The fact that the wave function is multi-valued means that you can't just apply the simple argument made in some textbooks that the phase acquired on exchange must square to 1.
Particle exchange corresponds to moving along closed loops in configuration space. If the loop is non-contractible, the wave function is free to pick up an arbitrary complex phase when the system moves along this curve. In 2D, the configuration space of two particles can be thought of as a cone with the apex point removed, so a path that circles the apex $n$ times is not deformable into one that circles it $m\neq n$ times. That is, a double exchange is not deformable into the trivial path. But in 3D, the configuration space of two particles (for fixed separation between the particles) is a sphere with antipodal points identified. This space is such that travelling twice around a non-contractible loop results in the trivial loop. That is to say, a double exchange is equivalent to no exchange. Hence any phase acquired by a single exchange must square to 1.
Further Details
Let's assume we have a wave function with no spin components. That is, at each point in configuration space, $\mathbf{x}\in M$, we have a vector $\mathbf{\Psi}(\mathbf{x})$ in a 1D complex Hilbert space, $h_\mathbf{x}= \mathbb{C}$. Written in terms of a normalised basis, $\mathbf{\chi}_\mathbf{x}$, $\mathbf{\Psi}(\mathbf{x}) = \psi (\mathbf{x}) \mathbf{\chi}_\mathbf{x}$. Seeing as the space is 1D, $\mathbf{\chi}_\mathbf{x}$ is just a complex phase, $e^{i\theta_\mathbf{x}}$. The function $\psi (\mathbf{x})$ is the wave function of the system and is multi-valued, while $\mathbf{\Psi}(\mathbf{x})$ is the state of the system in Hilbert space, and is assumed to be single valued. In other words, we have a fiber bundle where the configuration space is the base manifold and the fiber is $h_\mathbf{x}$. The state $\mathbf{\Psi}$ is a section of the fiber bundle.
Given a choice of basis $\mathbf{\chi}_\mathbf{x}$ at point $\mathbf{x}$, we can define a basis at all points in $M$ if we have a notion of parallel transport. Given a certain choice of (flat) connection, one finds no change to the basis after parallel transport along any curve, unless that curve circles around one of the punctures in the space. Along one of these non-trivial curves, $\mathcal{C}$, the basis vectors acquire a complex phase, $e^{i\alpha(\mathcal{C})}$, independent of $\mathbf{x}$. In order for the state $\mathbf{\Psi}$ to be single-valued, the wave function must be changed by the inverse complex phase to that acquired by the basis vector: $e^{-i\alpha(\mathcal{C})}$.
If any two paths are homotopic (deformable into each other), they must circle around these punctures the same number of times and so they acquire the same factor. That is, the phase acquired by exchange is governed by the fundamental group of the configuration space. Indeed, $e^{i\alpha(\mathcal{C})}$ belongs to a unitary representation of the fundamental group acting on the vector space $h_\mathbf{x}$.
In 3D, the fundamental group of configuration space is isomorphic to the symmetric group, since a double exchange is homotopic to the trivial element (no exchange). The only 1D representations of the symmetric group are the trivial rep (bosons) and the alternating rep (fermions), which sends odd permutations to -1 and even permutations to +1. Meanwhile in 2D, the fundamental group for 2 particles is isomorphic to $\mathbb{Z}$ since exchanging the particles $n$ times is distinct from doing so $m\neq n$ times. Thus the fundamental group for $N$ particles is isomorphic to the braid group. The braid group has a continuum of 1D representations corresponding to multiplication by any complex phase.
A: You should try to check that the $n$ particles lagrangian
$$
L = L_{0} - \alpha \sum_{i =1, j > i}^{n}\frac{d}{dt}\theta(\mathbf r_{i} - \mathbf r_{j}),
$$
where $L_{0}$ is "standard" many particles lagrangian and $\theta (\mathbf r_{i} - \mathbf r_{j})$ is azimuthal angle between $i$th and $j$th particles position vectors, describes the statistics of anyons (the amplitude given of particles intercharging acquires additional phase proportional to $\alpha$). 
Next, try to compute the Hamiltonian and ensure that anyons can be described by adding fictious non-propagating gauge field with some constrain.
Finally, the only thing which is remained to do is to solve corresponding Schroedinger equation...
P.S. About your comment ("...That is, we can only have fermions or bosons. We cannot have anyons..."): you don't take into account that the first homotopy group for 2-dimensional space is not the permutation group $Z_{2}$, as for spaces with dimension $> 3$. Therefore a loop encircled around the point in configuration space of particles living in 2D can't be shrinked to a point (unlike the 3D case). It's very simple to understand, if You imagine the closed loop with a point inside it in 2D and in 3D. This means, that the anyon's wave function need not be single valued under particles double interchange, and hence $\alpha $ from your comment need not be modulus 1.
A: To be more explicit, for future readers, there are second quantized operators that describe anyon 'excitation' or 'creation'. In the Kitaev model for example, he arrives at an effective Hamiltonian
\begin{equation}
H=-J_{e} \sum_{\text {vertices }} A_{s}-J_{m} \sum_{\text {plaquettes }} B_{p}, \quad \text { where } \quad A_{s}=\prod_{\text {star }(s)} \sigma_{j}^{x}, \quad B_{p}=\prod_{\text {boundary }(p)} \sigma_{j}^{z}.
\end{equation}
These $A,B$ operators give elementary $\mathbb{Z}_2$ $\mathbb{Z}_2$-excitations that give charges with energy $2J_e$ and vortices with charge $2J_m$. Those charges and vertices are exactly the $e,m$ Abelian anyons.
There is rarely a first quantized description of these particles for the same reason there is rarely a first quantized description of electrons in condensed matter models; local degrees of freedom such as absolute position are rarely important to describe dynamics.
A: Consider the wavefunction
\begin{equation}\label{eq:anyon_example}
 \Psi(z_1,z_2,\ldots, z_N)=\prod_{1\leq i<j\leq N} (z_i-z_j)^{p/q} e^{-\sum^N_{j=1} |z_j|^2},
\end{equation}
where we use complex coordinates $z=x+iy$ to label particle positions, and $p,q$ are relatively coprime constant integers. When $q\neq 1$ the wavefunction is multivalued, and we need to distinguish between clockwise (CW) and counter-clockwise (CCW) exchanges. For any pair $i,j$ of particles, if we braid  $z_i$ CCW  around $z_j$, we have $\Psi\to e^{i\pi p/q}\Psi$ while for a CW braid, we have $\Psi\to e^{-i\pi p/q}\Psi$, so this represent a wavefunction of $N$ identical abelian anyons.
