Conventional fMRI data analysis uses a model of the hemodynamic response based on a gamma function. The form of the fMRI response to stimuli can be modeled by convolution of the hemodynamic response, which related to the presented stimulus. This approach is based on the linear system analysis. However, it is known that most biological system should be nonlinear. It is also known that the hemodynamic response is varied due to subjects and brain regions. In this work, we hypothesized that all the neural activity of the human cerebral region measured by fMRI is not underlying linear system-based hemodynamic response. We investigated the hemodynamic response of the fusiform face area (FFA) and the primary visual area (V1) using the estimated $1^{st}$ and $2^{nd}$ order Volterra kernels based on system identification using Laguerre expansion technique. The results of analysis using the estimated Volterra kernels from block-designed fMRI data showed significantly more activation compared to the results of analysis using the hemodynamic response modeled by a gamma function. Furthermore, we could show the different characteristic of the hemodynamic responses between primary visual area and the fusiform face area, and the different characteristics of the hemodynamic responses due to subjects with the estimated $1^{st}$ and $2^{nd}$ order Volterra kernel based on system identification using Laguerre expansion technique.