We report detection of PSR B0656+14 with the Gran Telescopio Canarias in narrow optical F657, F754, F802, and F902 and near-infrared JHKs bands. The pulsar detection in the Ks band extends its spectrum to 2.2 $\mu$ m and confirms its flux increase towards the infrared. We also present a thorough analysis of the optical spectrum obtained by us with the VLT. For a consistency check, we revised the pulsar near-infrared and narrow-band photometry obtained with the HST. We find no narrow spectral lines in the optical spectrum. We compile available near-infrared-optical-UV and archival 0.3-20 keV X-ray data and perform a self-consistent analysis of the rotation phase-integrated spectrum of the pulsar using unified spectral models. The spectrum is best fitted by the four-component model including two blackbodies, describing the thermal emission from the neutron star surface and its hot polar cap, the broken power law, originating from the pulsar magnetosphere, and an absorption line near ∼0.5 keV detected previously. The fit provides better constraints on the model parameters than using only a single spectral domain. The derived surface temperature is $T_{NS}^{\infty } = 7.9(3)\times 10^5$ K. The intrinsic radius (7.8-9.9 km) of the emitting region is smaller than a typical neutron star radius (13 km) and suggests a non-uniform temperature distribution over the star surface. In contrast, the derived radius of the hot polar cap is about twice as large as the 'canonical' one. The spectrum of the non-thermal emission steepens from the optical to X-rays and has a break near 0.1 keV. The X-ray data suggest the presence of another absorption line near 0.3 keV.