Perl 5 + -pF//?(dd)/ -Mbigint, 74 bytes

$_="1357
7531
3175
5713";s!.!$F[$&]-3+(3333420615514242=~/./g)[$i++].$"!ge

Try it online!

Verify all test cases online!

Explanation

The -F flag, splits the input into @F as ('',22,'',12,'',18,'',87). In the code $_ is set to a string containing the index of @F that we need the number from, then s/// replaces any char (. - this excludes newline) with the the corresponding value from @F minus 3, plus the digit from 3333420615514242 at the $index, appending $" (space).

perl -nl -M5.010, 106 bytes

/.(dd).(dd)/;$,=$";say$`,$1,$2,$';say$'+1,$2-1,$1-3,$`+3;say$1-2,$`+2,$'+2,$2-2;say$2+1,$'-1,$`+1,$1-1

Try it online!

How does it work?

/.(dd).(dd)/;

This lines parses the input. It matches the /MM/yy part of the input, leaving MM in $1 (the first capture) and yy in $2. DD is available in $` (the "prematch"), and YY is available in $' ("postmatch").

$, = $";

This sets the output field separator equal to the list separator, the latter is by default a space. This makes that our output numbers are separated by a space.

say $`    , $1    , $2    , $';
say $' + 1, $2 - 1, $1 - 3, $` + 3;
say $1 - 2, $` + 2, $' + 2, $2 - 2;
say $2 + 1, $' - 1, $` + 1, $1 - 1

Print out the 16 required numbers. In the top row, it will keep any leading spaces, as the values are printed "as is".

JavaScript (Node.js),  88  84 bytes

Saved 4 bytes thanks to @NahuelFouilleul

d=>`-BWe
fV?0
@/gU
Xd.A`.replace(/./g,c=>d.substr((n=Buffer(c)[0])/7-6,2)-3+n%7+' ')

Try it online!

How?

Each value in the matrix is encoded as a character whose ASCII code is:

$$7k+o+42$$

where $k$ is the 0-indexed position of the variable in the input string and $o-3$ is the offset to add.

Example:

For MM-2, the position is $3$:

"DD/MM/yyYY"
    ^
 0123456789

and the offset is$-2$, leading to:

$$7times3+(-2+3)+42=64$$

which is '@'.

Jelly,  27  26 bytes

ḟṂs2Vṙ1+J$_Ɱ4F“Œñḳ⁶.C’œ?ṁƲ

A monadic Link accepting the list of characters as specified which yields a list of lists of integers.

Try it online!

How?

ḟṂs2Vṙ1+J$_Ɱ4F“Œñḳ⁶.C’œ?ṁƲ - Link: list of characters, D    e.g. "01/09/2000"
 Ṃ                         - minimum                             '/'
ḟ                          - filter discard                      "01092000"
  s2                       - split into twos                     ["01","09","20","00"]
    V                      - evaluate                            [1,9,20,0]
     ṙ1                    - rotate left 1                       [9,20,0,1]
         $                 - last two links as a monad:
        J                  -   range of length                   [1,2,3,4]
       +                   -   add                               [10,22,3,5]
           Ɱ4              - map across [1..4] with:
          _                -   subtract                          [[9,21,2,4],[8,20,1,3],[7,19,0,2],[6,18,-1,1]]
                         Ʋ - last four links as a monad - f(x):
             F             -   flatten                           [9,21,2,4,8,20,1,3,7,19,0,2,6,18,-1,1]
              “Œñḳ⁶.C’     -   base 250 integer = 19644039699318
                      œ?   -   nth permutation                   [1,9,20,0,1,19,6,4,7,3,2,18,21,-1,2,8]
                        ṁ  -   mould like (x)                    [[1,9,20,0],[1,19,6,4],[7,3,2,18],[21,-1,2,8]]

Charcoal, 46 bytes

Ｅ⪪+5?IJ>2.3-K=@H,4⁴⪫Ｅι◧Ｉ⁺⁻﹪℅λχ³§Ｉ⪪⪫⪪θ/ω²÷℅λχ² 

Try it online! Link is to verbose version of code. Explanation:

Ｅ⪪+5?IJ>2.3-K=@H,4⁴

Split the string literal +5?IJ>2.3-K=@H,4 into chunks of 4 characters and loop over each chunk, which represents a row of the magic square.

⪫Ｅι◧Ｉ...² 

Loop over each character, and pad each result to 2 characters and join them together to form the row.

⁺⁻﹪℅λχ³§...÷℅λχ

Take the ASCII code of the character and divmod it by 10. Subtract 3 from the remainder to give a value in the range -3..3 and add the date part given by the quotient (this is offset by 4 but that works due to Charcoal's cyclic indexing).

Ｉ⪪⪫⪪θ/ω²

Delete the /s from the input string, split the result into pairs of digits and take their value.

K (ngn/k), 53 51 bytes

Solution:

{4 4#x,(3-75330941200)+(x:.'4 2#x^"/")@414962353}

Try it online!

Explanation:

Probably further golfable...

{4 4#x,(3-75330941200)+(x:.'4 2#x^"/")@414962353} / the solution
{                                                 } / lambda taking implicit x
                                        414962353  / decode 14962353 into base 4 => 3 2 1 0 1 0 3 2 2 3 0 1
                                       @            / use that to index into
                        (             )             / do all this together
                                 x^"/"              / x except (^) "/"
                             4 2#                   / reshape as 4x2 grid
                           .'                       / 'value each', convert e.g. "20" to 20
                         x:                         / save back into x
                       +                            / add to
       (              )                             / do this together
          75330941200                              / decode 5330941200 into base 7 => 2 4 6 0 5 1 1 5 2 4 2 4
        3-                                          / subtract from 3 => 1 -1 -3 3 -2 2 2 -2 1 -1 1 -1
     x,                                             / prepend x
 4 4#                                               / reshape into 4x4 grid

Extra:

• -1 byte with ^"/" (except "/") instead of dropping indices 5 and 2 (_/5 2)
• -1 byte with 3-75330941200 instead of -3+78510346000

APL (Dyalog Unicode), 66 bytes

{4 4⍴(∊(⊣,-)¨3-⍨1011304⊤⍨8⍴6)+,(⊢⍪⌽∘⊖⍪⌽⍪⊖)2 2⍴10⊥⍉⍎¨4 2⍴'/'⎕R''⊢⍵}

Try it online!

How?

The key insight is that the rows correspond to a single 2×2 matrix that has been mirrored vertically, horizontally, not at all, or both.

{4 4⍴(∊(⊣,-)¨3-⍨1011304⊤⍨8⍴6)+,(⊢⍪⌽∘⊖⍪⌽⍪⊖)2 2⍴10⊥⍉⍎¨4 2⍴'/'⎕R''⊢⍵} ⍝ Anonymous inline function
  ⍵            ⍝ Input, e.g. '22/12/1887'
  '/'⎕R''      ⍝ Replace '/' with '', giving e.g 22121887
  10⊥⍉⍎¨4 2⍴   ⍝ Reshape and execute to give vector (22,12,18,87)
  2 2⍴         ⍝ Reshape to 2×2 matrix (22 12)(18 87)
  (⊢⍪⌽∘⊖⍪⌽⍪⊖)  ⍝ Mirror the matrix in 4 ways to get the 4 rows
  ,            ⍝ Flatten to a row vector
  1011304⊤⍨8⍴6 ⍝ Base 6 encoding of 3 3 4 0 1 5 4 4
  3-⍨          ⍝ Subtract 3 from each: 0 0 1 ¯3 ¯2 2 1 1
  ∊(⊣,-)¨      ⍝ Append the negative to each and flatten: 0 0 0 0 1 ¯1 ¯3 3 ¯2 2 2 ¯2 1 ¯1 1 ¯1
  +            ⍝ Add these constants to each element
  4 4⍴         ⍝ Reshape to give a 4×4 matrix  

05AB1E, 27 26 24 bytes

þ2ôÀā+āδ-˜•Itíˀߕ.I4ô

-2 bytes by porting @JonathanAllan's Jelly answer, so make sure to upvote him as well!

Try it online or verify all test cases.

Old 27 26 bytes approach:

þ2ôÂD2ôR˜Â)•CV¤½‘àý•S4ä3-+

Try it online or verify all test cases.

Explanation:

þ         # Leave the digits of the (implicit) input-string: "dd/MM/YYyy" to "ddMMYYyy"
 2ô       # Split it into parts of size 2: [DD,MM,yy,YY]
À         # Rotate it once towards the left: [MM,yy,YY,DD]
 ā        # Push a list in the range [1,length] (without popping): [1,2,3,4]
  +       # Add the values at the same indices: [MM+1,yy+2,YY+3,DD+4]
ā         # Push [1,2,3,4] again
 δ        # Apply double-vectorized on these lists:
  -       #  Subtract
          # (results in: [[MM,MM-1,MM-2,MM-3],[yy+1,yy,yy-1,yy-2],[YY+2,YY+1,YY,YY-1],
          #  [DD+3,DD+2,DD+1,DD]])
   ˜      # Flatten it:
          #  [MM,MM-1,MM-2,MM-3,yy+1,yy,yy-1,yy-2,YY+2,YY+1,YY,YY-1,DD+3,DD+2,DD+1,DD]
    ā     # Push a list in the range [1,length] (without popping) again: [1,...,16]
•IZΔmλR•  # Push compressed integer 19556972234337
 .I       # Get the 19556972234337'th permutation of the [1,16] list
   è      # Index each into the list of values (modular 0-based, so 16 is the first):
          #  [DD,MM,yy,YY,YY+1,yy-1,MM+3,DD+3,MM-2,DD+2,YY+2,yy-2,yy+1,YY-1,DD+1,MM-1]
    4ô    # And split it into parts of size 4 again
          # (after which the result is output implicitly

þ2ô       # Same as above: [DD,MM,yy,YY]
         # Bifurcate this list (short for Duplicate & Reverse copy)
          #  [DD,MM,yy,YY] → [YY,yy,MM,DD]
 D        # Duplicate it
  2ô      # Split it into pairs
          #  [YY,yy,MM,DD] → [[YY,yy],[MM,DD]]
    R     # Reverse the pairs
          #  [[YY,yy],[MM,DD]] → [[MM,DD],[YY,yy]]
     ˜    # Flatten it back again to a quartet
          #  [[MM,DD],[YY,yy]] → [MM,DD,YY,yy]
      Â   # Bifurcate that as well
          #  [MM,DD,YY,yy] → [yy,YY,DD,MM]
       )  # And wrap all list on the stack into a list
          #  → [[DD,MM,yy,YY],[YY,yy,MM,DD],[MM,DD,YY,yy],[yy,YY,DD,MM]]
•CV¤½‘àý• # Push compressed integer 3333420615514242
 S        # Convert it to a list of digits: [3,3,3,3,4,2,0,6,1,5,5,1,4,2,4,2]
  4ä      # Split it into 4 equal-sized parts: [[3,3,3,3],[4,2,0,6],[1,5,5,1],[4,2,4,2]]
    3-    # Subtract 3 from each: [[0,0,0,0],[1,-1,-3,3],[-2,2,2,-2],[1,-1,1,-1]]
      +   # Add the values in the lists together:
          #  [[DD,MM,yy,YY],
          #   [YY+1,yy-1,MM-3,DD+3],
          #   [MM-2,DD+2,YY+2,yy-2],
          #   [yy+1,YY-1,DD+1,MM-1]]
          # (after which the result is output implicitly

See this 05AB1E tip of mine (sections How to compress large integers? and How to compress integer lists?) to understand why •IZΔmλR• is 19556972234337 and •CV¤½‘àý• is 3333420615514242.