random walk, how many steps does it take to visit all the states? For example, random walk in a 2-D integer space, 
range of x: [-5, 5], range of y: [-8, 8]
starting at (0,0), 
equal probability to go left/right/up/down, step size 1.
How many steps it is going to take to visit all the states?
Is there a math equation to calculate this statistically?
 A: This is not really an answer, more a thought on how to get started, but it's too long for a comment.  For small rectangle, it's possible to solve the problem exactly, though the method I have in mind quickly becomes unwieldy.  If we have a $w\times h$ rectangle, we can model the problem as a finite-state absorbing Markov chain with states of the form $(j,k,s)$ where $0\leq j\leq w,\ 0\leq k\leq h,$ and $s$ is the set of states not yet visited.  We can compress all the states with $s = \emptyset$ into a single, absorbing state.  Then the value you want is the expected time to absorption.
You might try this for small values of $w,h$ and see what you get.  At the number of states is $wh(2^{wh-1}+1),$ this will only work for very small values of $w$ and $h$.  
I suggest starting by trying to get an answer for $h=1.$  A big advantage is that the problem is one-dimensional, so the we can represent the states as $(x, l, r)$ where $x$ is the current position, and $l$ and $r$ are the furthest left and right that we have traveled, respectively.  This makes the number of states much more manageable.
I wrote a python script to compute the expected times in the one-dimensional case.  You don't specify the behavior when the walk reaches the boundaries.  I assumed that it was at the left (respectively right) boundary, it would stay put with probability $1/2$ and move right (respectively left) with probability $1/2.$  I only print the expected time when the walk starts in the center of the range, but the script actually computes the expected time from all states.
Here is the script:
import numpy as np

def makeTransitions(R):
    states = [(pos,left,right) for left in range(R+1) 
                       for right in range(left, R+1) 
                       for pos in range(left,right+1) 
                       if left != 0  or right != R]
    table = { }
    for idx, s in enumerate(states):
        table[s] = idx
    n = len(states)
    P = np.zeros((n,n))
    for s in table:
        pos,left,right = s
        s1 = min(pos+1,R), left, min(max(pos+1,right), R)
        s2 = max(pos-1,0), max(0, min(pos-1,left)), right
        for t in (s1, s2):
            if (t[1], t[2]) == (0, R):  
                continue    # absorbed
            P[table[s], table[t]] = .5
    return table, P

def time(R):
    states, P = makeTransitions(R)
    n = len(states)
    I = np.eye(n)
    N = np.linalg.inv(I-P)
    J = np.ones((n,1))
    E = N@J
    start = states[(R//2,R//2,R//2)]
    return E[start][0]

for R in range(1,15):
    print(R, time(R))

This produces the output:
1 2.0
2 7.0
3 13.999999999999998
4 23.999999999999993
5 35.99999999999999
6 50.99999999999998
7 68.0
8 87.99999999999999
9 110.0
10 135.0
11 162.0
12 192.00000000000003
13 224.0
14 259.0

Allowing for rounding error, this agrees with OEIS A249547.  So far, I am at a loss to see how to prove that the pattern hold in general.
