What is difference between fermions and spins? A spin model i.e. $H_s = \sum_i^{L-1} S_i^x\cdot S_{i+1}^x$  can be written in matrix form as following
$$H_s = \big(S_1^x \otimes S_2^x \otimes I_3^2 \otimes I_4^2\otimes \cdots\otimes I_{L-1}^2\big) 
+ \big( I_1^2 \otimes S_2^x \otimes S_3^x \otimes I_4^2 \otimes  \cdots I_{L-1}^2\big) 
 + \cdots +\big( I_1^2 \otimes I_2^2 \otimes   I_3^2 \cdots \otimes S_{L-2}^x \otimes S_{L-1}^x   \big) 
$$ 
where $S_i^x =\begin{bmatrix} 0&0.5\\0.5&0\end{bmatrix} $ and $I_i^2 =\begin{bmatrix} 1&0\\0&1\end{bmatrix}$
On the other hand, let we have a fermionic model $H_f =  \sum_i^{L-1} c_i^\dagger c_{i+1}$. For spinless fermions, one can represent creation operator $c_i^\dagger$ and annihilation operator $c_i$ in matrix form as
$c_i^\dagger =\begin{bmatrix} 0&1\\0&0\end{bmatrix} $ and $c_i =\begin{bmatrix} 0&0\\1&0\end{bmatrix}$.
So if we can write fermionic operators in 2-by-2 matrix just like spin operators, we can write complete Hamiltonian matrix just like above equation i.e.
$$H_f = \big(c_1^\dagger \otimes c_2 \otimes I_3^2 \otimes I_4^2\otimes \cdots\otimes I_{L-1}^2\big) 
+ \big( I_1^2 \otimes c_2^\dagger \otimes c_3 \otimes I_4^2 \otimes  \cdots I_{L-1}^2\big) 
 + \cdots +\big( I_1^2 \otimes I_2^2 \otimes   I_3^2 \cdots \otimes c_{L-2}^\dagger \otimes c_{L-1}   \big) 
$$ 
Question:
If my approach is correct then what is the difference between spins models and fermionic models? Spin models obey commutation relations $[S_i^a,S_j^b] = i\hbar \epsilon_{abc} S_c \delta_{ij}$ while fermions obey anti-commutation relations $\{c_i^\dagger, c_j\} = \delta_{ij}$. Where do these relations go when we write Hamiltonian in matrix form?
 A: Your matrix expressions for the fermionic operators are wrong because they do not obey the anti-commutation relations. More precisely, they are correct if you have only a single fermionic mode, but are wrong for $L>1$. If you want to get a matrix representation of fermionic operators you need to use the Jordan-Wigner transformation:
\begin{align}\hat{c}_k & =  \hat{\sigma}_k^- \prod_{j<k}\hat{\sigma}^z_j \\
& = \bigotimes_{j=1}^{k-1}\hat{\sigma}^z_j \otimes \hat{\sigma}^-_k \bigotimes_{l = k+1}^L \hat{I}^2_l, 
\end{align}
where the second line is the tedious full expression making the tensor product structure explicit. Here, $\hat{\sigma}_k^- = (\hat{\sigma}_k^x - {\rm i}\hat{\sigma}^y_k)/2$ is the spin lowering operator on site $k$. From this you can construct any fermionic operator you want as strings of $\hat{c}_k$, $\hat{c}_l^\dagger$, and linear combinations thereof.
As a simple example of how this works, consider the simplest non-trivial case of $L=2$ sites. If one (wrongly!) postulates that $\hat{c}_{1} = \hat{\sigma}^-_1\otimes\hat{I}^2_2$ and $\hat{c}_{2} = \hat{I}^2_1\otimes\hat{\sigma}^-_2$, the anti-commutation relation would be $\{\hat{c}_1,\hat{c}_2\} = 2\hat{\sigma}^-\otimes\hat{\sigma}^-\neq 0$, so that's no good. Using the correctly defined operators instead, we have
\begin{align}
\{\hat{c}_1,\hat{c}_2\} &= \hat{\sigma}^-_1 \hat{\sigma}^-_2 \hat{\sigma}^z_1 + \hat{\sigma}^-_2\hat{\sigma}^z_1\hat{\sigma}^-_1\\
& = \hat{\sigma}^- \hat{\sigma}^z \otimes \hat{\sigma}^- + \hat{\sigma}^z\hat{\sigma}^- \otimes \hat{\sigma}^- \\
& = \{\hat{\sigma}^-,\hat{\sigma}^z\} \otimes \hat{\sigma}^- \\
& = 0,
\end{align}
since $\{\hat{\sigma}^-,\hat{\sigma}^z\} = 0$.
