Why do we need large system in Monte Carlo simulations? Why do we need to simulate large systems (especially in spin systems)? Is it just a better representation of the real behaviur of magnetic systems? Is there an instability of the finite system with respect to phase transitions?
Finally, what does it mean to reach the thermodynamic limit for a spin glass system?
 A: As @zivo stated in the comment, the primary motivation is to avoid "edge" effects. If you arrange $1000$ spins in a $10\times10\times10$ cube, then $488$ of them (nearly half) lie on the faces of the cube, and $512$ in the interior. If we go to $100\times100\times100=10^6$ spins, the fraction on the cube surface is still $\sim 6$%.
However that's not the end of the story. Almost all simulations of bulk systems use periodic boundary conditions which eliminate the physical surfaces and make the system translationally invariant. So every spin is equivalent to every other spin. However, this does not completely get rid of finite size effects. The system is artificially periodic: if the side of the cubic box is $L$, then the system is effectively periodic with period $L$ in all three directions. Another way of looking at this is, if you are interested in describing properties of the system in Fourier-transformed variables, the components of the wave-vector $\mathbf{k}$ are restricted to integer multiples of $k_{\text{min}}=2\pi/L$. No behaviour having a longer wavelength than $L$ can be modelled. In the real system, there is no lower limit to the wave-vector components, and no upper limit to the wavelengths of physical phenomena.
This may not be important for many applications of Monte Carlo simulation, but is likely to be most important near phase transitions, where correlation lengths and fluctuation quantities such as susceptibilities and heat capacities diverge. However, it turns out that in most cases the effects of finite size on such properties are quite well understood, in the vicinity of first-order phase transitions, and continuous ones (where the critical exponents for various universality classes are known). 
I would not describe it as an "instability", as in your question.
Instead, for first-order phase transitions,
the discontinuous or sharp changes that you expect in certain thermodynamic quantities get "smoothed out" and "broadened" by the finite-size effects,
and the transition temperature (say) gets shifted a bit.
Those effects are typically of order $1/N$ where $N$ is the number of spins.
For continuous transitions, 
the divergent correlation lengths are limited by the finite size,
and this affects the various scaling laws in a predictable way.
This is the area of finite size scaling theory, and there are a few general descriptions of it in statistical mechanics text books (for example Scaling and renormalization in statistical physics by John Cardy, Cambridge University Press, 1996) and online.
So a typical scenario is to study several system sizes, and scale the results using the anticipated powers of $L$, to make all the plots of simulation results collapse onto a unified curve. This kind of approach was pioneered by Kurt Binder (see this page) in the 1980's, and he has many many publications using this approach. It means that we don't need to simulate an infinitely large system to understand the details of critical phenomena!
For spin glasses, things get complicated. Often, there is a random element to the microscopic specification of the model: for example, the coupling constants. One often adopts the idea of the "quenched" average: computing simulation averages from suitably long runs, for a given set of coupling constants. The thermodynamic limit then implies averaging over the possible random choices of coupling constants, afterwards. On top of this, the variation with system size should be studied, and the large-$L$ behaviour of a spin glass may be nontrivial. 
So, I can't give a definitive answer to this part of your question, but I can refer you to these lecture notes by Peter Young (which also explains how finite-size scaling can be used), and to the first two references therein: K Binder and AP Young, Rev Mod Phys, 58, 801 (1986) and AP Young (ed) Spin glasses and random fields (World Scientific, 1998).
