The Monte Carlo method was used to generate canonical-ensemble isotherms for N=25, 64, 100, and 400 parallel hard squares in a two-dimensional periodic box. The Monte Carlo realizations with their chain lengths ranging up to ∼ 8× 105 trials per particle were employed to ascertain statistical accuracy of the calculated pressures. The resulting isotherms are monotone increasing functions of the density with a possibility of a higher-order transition but without showing any sign of a first-order phase transition, with a density increment greater than 0.005ρ0 (ρ0 is the close-packed density). This departs remarkably from the behavior of the hard-disk isotherm, which displayed a van der Waals-like loop for N ≳ 72. Instead, it resembles isotherms of the lattice-gas systems (e.g., a square-lattice system with the first and second nearest-neighbor exclusion, or a system of dimers) which, like the present case, contain a residual degree of freedom at close packing. By numerically integrating the Monte Carlo isotherms, the constant which corrects the entropy expression of the self-consistent free-volume theory at close packing was evaluated. The resulting value (2 ln2+0.268± 0.004) for the N=400 case agrees satisfactorily with the corresponding value (2 ln2+0.26042) predicted from the cell-cluster theory for the checkerboard structure. The N dependence of the Monte Carlo isotherms was examined in both the low- and the high-density regions. In the low-density region, the Monte Carlo data are in satisfactory agreement with theoretically predicted values, whereas, in the high-density region, the MC data revealed existence of the N -dependent contribution to the pressure attributable to the correlated slidings of rows or columns of particles. This can occur even at 95% of the close-packed density for N≤ 100.