For example, if I am using neural network for regression analysis, neural network without any hidden layer will give linear function. Does hidden layer corresponds to degree of approximation function? If so, how neural network with single hidden layer can be a universal approximator ?
It can be true that more layers/depth can yield more nonlinearity, but I do not think using the analogy of a polynomial's degree is a perfect analogy and it's certainly not mathematical fact.
You can look at deep net regression a few different ways: learning a nonlinear function. In this sense, it is very natural to think in terms of degrees and polynomials. Another sense is that learning using neural networks involves learning a partition of the domain into sub-domains, and on each sub-domain we approximate the function with a linear superposition or a deep composition of simpler functions. In this sense, we measure representational capacity in terms of how many sub-domains we can learn. Of course, these concepts may not be orthogonal, but they emphasize different aspects. The polynomial degree analogy feels a lot like interpolation to me while I believe the philosophy (and recent theory) with deep net representational capacity is more about partitioning or learning regions.
There's a couple of different angles to your question. If you are asking why is a single hidden layer NN is a Universal Function Approximator, that's a mathematical question, and the proof is available, for example see Cybenko's proof for sigmoidal nonlinearities. One key quote from Cybenko's paper is the following:
"that networks with one internal layer and an arbitrary continuous sigmoidal function can approximate continuous functions wih arbitrary precision providing that no constraints are placed on the number of nodes or the size of the weights"
So to approximate a continuous function $f$ with a single hidden layer NN to some given precision $\epsilon$ may require an extremely wide architecture, despite being incredibly shallow. Meaning, we have many many many different sigmoidal activation functions. By allowing ourselves to have possibly many such sigmoidal activation functions, i.e. the wide architecture, we are allowing ourselves to have a much more intricate partition of our domain, which increases our approximation accuracy.
There's some further ways to break down the Universal approximation theorem. For example, it requires the activation function be nonconstant, bounded, and monotonically-increasing continuous function. Take the case of the sigmoidal activation, which asymptotes out to 0 as $x \rightarrow -\infty$ and to 1 as $x \rightarrow +\infty$. One of the original motivations for using the sigmoidal activation was that it can act like a switch. Allowing the weights to have arbitrary magnitude, i.e. no hard limit on the size of the parameter values, we can effectively shift sigmoids into their "off" states with very negative values and shift sigmoids into their "on" states with very positive values. And training a single hidden layer corresponds to learning a good configuration of parameters. By allowing for the weights to have unbounded size (both in the negative and positive sense), we can interpret the single hidden layer NN as partitioning the domain into sub-spaces where a specific configuration of the sigmoidals are "on" and contribute to the function approximation and the others are switched "off." Now if we allow ourselves to have a ton of these sigmoids, you start to get some intuition as to why the sigmoidal Universal Approximation Theorem is true.
One question that I believe has motivated researchers is: how efficient is this? In particular: how does representational capacity scale with the number of hidden units for single hidden layer neural networks? For deeper architectures?
I believe recent work by Yoshua Bengio and his collaborators have taken these ideas further with the wildly successful ReLU $f(x) = \max(0, x)$ and a generalization called the Maxout activation functions. These ReLU and Maxout networks tend to have multiple layers and are not shallow. Please see the following How to understand the geometric intuition of the inner workings of neural networks? for a more detailed explanation. But roughly speaking, they proved that Maxout networks learn an exponential number of linear regions. In particular, the number of linear regions learned scales exponentially with the number of layers.
So perhaps thinking about degree in terms of polynomials is not the best analogy, but it's better to think about depth in terms of how it scales with the number of regions (size of the partition) our network can learn.
In contrast to single hidden layer NN, if the general principle is that representational capacity scales with depth, then we would require an extremely wide and shallow architecture and on a per-parameter basis, this may be very inefficient compared to a deeper network.