Multiply Pauli Matrices
Python 2, 108 89 87 86 bytes

x=y=0
for m in map(int,raw_input()):x+=m*y and(m-y)%3*3/2;y^=m
print"--i"[~x%4::2]+`y`


(Thanks to @grc and @xnor for the help)

Explanation

Let's split up the coefficient and the base matrix. If we focus on the base matrix only, we get this multiplication table (e.g. 13 is -i2, so we put 2):

  0123

0 0123
1 1032
2 2301
3 3210


which just happens to be the same thing as doing bitwise xor.

Now let's focus on the coefficients. If we let 0123 denote 1,i,-1,-i respectively, we get:

  0123

0 0000
1 0031
2 0103
3 0310


For this we first check if either number is 0 by doing m*y, taking care of the left column and top row. Adding in (m-y)%3 then gives:

  0123

0 0000
1 0021
2 0102
3 0210


which is close, except that we have 2 instead of 3. We account for this by performing *3/2.

For indexing, we notice that if we take the string "--i" and select every second character starting from indices 0123 we get "-i","-","i","".